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ABSTRACT 


Richmond's  moment-method  analysis  for  the  current  distribution  and  input  impedance  of 
a  monopole  element  on  a  disk  ground  plane  above  flat  Earth  is  used  to  obtain  the  far-zone 
field  in  the  free-space  region.  Numerical  results  for  directivity  and  radiation  efficiency  are 
presented  as  separate  entities,  unlike  previously  reported  results  based  on  Monteath's 
compensation  theorem  or  Sommerfeld's  attenuation  function  that  give  only  the  product  of  the 
directivity  and  radiation  efficiency. 
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SECTION  1 


INTRODUCTION 


The  modeling  of  monopole  elements  with  circular  ground  planes  in  proximity  to  Earth 
has  been  greatly  enhanced  in  recent  years  by  method -of-moments  programs  developed  by 
Richmond  for  disk  ground  planes  [1]  and  by  Burke,  et  al„  for  radial  wire  ground  planes 
[2,3,4]. 

The  method-of-moments  models  include  the  following  advantages  over  models  based  on 
Monteath's  compensation  theorem  [5,6,7,8,19]  or  Sommerfeld's  attenuation  function  [9]: 

(1)  current  on  the  ground  plane  is  computed  rather  than  approximated  by  that  for  a  perfect 
ground  plane;  (2)  results  are  valid  not  only  for  moderately  large  ground  planes,  but  for 
electrically  small  ground  planes;  (3)  ground-plane  edge  diffraction  is  determined  directly 
rather  than  neglected  or  obtained  by  perturbation  methods;  (4)  analytical  restrictions  on 
evaluating  Sommerfeld's  integral  (such  as  requiring  that  the  Earth's  complex  relative 
permittivity  have  a  modulus  much  greater  than  unity)  are  avoided;  and  (5)  directivity  and 
radiation  efficiency  are  determined  as  separate  entities,  rather  than  being  lumped  together  as 
a  product  to  yield  the  antenna  gain.  Nevertheless,  those  other  models  are  useful  for 
validating  method-of-moments  numerical  results  and  for  treating  large  ground  planes  whose 
segmentation  in  method-of-moments  models  would  exceed  computer  computational  capacity 
and  precision. 

Richmond  has  presented  a  moment-method  analysis  for  the  current  distribution  and  input 
impedance  of  a  monopole  element  on  a  disk  ground  plane  in  free-space  [10]  and  above  flat 
Earth  [1]  with  numerical  evaluation  by  computer  programs  RICHMD1  and  RICHMD3, 
respectively.  Weiner,  et  al.  [1 1,12]  have  used  Richmond's  results  in  reference  10  to  develop 
a  computer  program,  RICHMD2,  for  the  far-zone  field,  directivity,  and  radiation  efficiency 
for  the  case  when  the  ground  plane  is  in  the  free  space.  The  present  effort  uses  Richmond's 
results  in  reference  1  to  obtain  the  far-zone  field  when  the  ground  plane  is  above  flat  Earth. 
Numerical  evaluation  of  the  far  zone  field  is  achieved  with  Richmond's  computer  program 
RICHMOND4. 
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Consideration  is  limited  to  the  far- zone  field  in  the  free-space  region,  with  ionospheric 
effects  excluded.  When  the  observer  approaches  the  air-Earth  interface,  the  total  far- zone 
field  will  include  a  small  contribution  from  the  "surface  wave,"  but  this  term  is  not 
considered  here.  Instead  of  a  null  on  the  radio  horizon,  the  surface  wave  will  contribute  a 
far- zone  evanescent  field  that  is  small  compared  to  the  far- zone  field  in  the  direction  of  peak 
directivity.  However,  the  significant  effect  of  the  near-field  surface  wave  in  reducing 
radiation  efficiency  is  included  in  the  present  analysis. 

The  radiation  field  from  the  surface  magnetic  current  density  (magnetic  frill)  of  the 
coaxial  line  feed  is  included  in  the  present  analysis.  Although  the  magnetic  frill  is  the 
excitation  source  for  the  current  on  the  monopole  element  and  the  disk  ground  plane,  its 
contribution  to  the  far- zone  field  may  usually  be  neglected  (as  was  done  in  references  1 1  and 
12)  unless  the  monopole  element  is  so  short  that  the  radiation  resistance  of  the  magnetic  frill 
becomes  comparable  to  that  of  the  monopole  element. 

The  theoretical  model,  numerical  results,  validation  of  numerical  results,  and  conclusions 
are  given  in  sections  2,  3, 4,  and  5,  respectively.  A  more  comprehensive  treatment  and 
review  of  monopole  elements  on  circular  ground  planes  in  proximity  to  flat  Earth  is  given  in 
reference  13. 
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SECTION  2 


THEORETICAL  MODEL 


2.1  METHODOLOGY 

The  antenna  geometry  consists  of  a  vertical  monopole  element  (length  h  and  radius  b),  on 
an  infinitely  thin  disk  ground  plane  of  radius  a  at  a  height  z0  above  flat  Earth  (see 
figure  1).  The  Earth,  with  a  dielectric  constant  £r,  conductivity  <7  (S/m)  at  a  radian  frequency 
to  (rad/s),  and  free-space  wavelength  X  (m),  has  a  complex  relative  permittivity  £*/£o 
=  £r  (1  -  j  tan  S )  where  tan  5  =  loss  tangent  =  o/((o  er  So)  =  {Xoftn  er )  ipo^^o)1^ 

«  60  Xa/er.  The  monopole  element  and  disk  are  assumed  to  have  infinite  conductivity.  The 
location  of  an  arbitrary  far-zone  observation  point  P  is  designated  by  spherical  coordinates 
ip,  6,  <P)  with  original  O  at  the  air-Earth  interface  below  the  monopole  element. 

The  feed  for  the  monopole  antenna  is  a  coaxial  line  with  its  inner  conductor  connected 
through  a  hole  of  radius  bj  in  the  ground  plane  to  the  vertical  monopole  element  and  its  outer 
conductor  connected  by  means  of  a  flange  to  the  ground  plane.  The  inner  conductor’s 
diameter  is  equal  to  the  monopole  element's  diameter  2b  and  the  outer  conductor's  diameter  is 
equal  to  the  ground-plane  hole  diameter  2bj.  The  current  on  the  outside  of  the  coaxial-line 
feed  is  assumed  to  be  zero  because  of  the  attenuation  by  lossy  ferrite  toroids  along  the 
exterior  of  the  coaxial-line  feed  (see  section  2.4  of  reference  12).  The  coaxial  line  feed 
excitation  may  be  replaced  by  an  equivalent  surface  magnetic  current  density  (magnetic  frill) 
Mq  given  by  equation  (23)  of  section  2.4. 

The  magnetic  frill  excitation  gives  rise  to  a  monopole  element  current  distribution  / Z(z) 
along  the  z  axis  of  the  monopole  element  and  a  disk  current  density  distribution  Jp{p )  in  the 

radial  direction  p  in  the  plane  of  the  disk.  The  current  density  Jp(p )  is  the  net  current 
density  on  the  top  and  bottom  of  the  disk  (see  equation  2.4.1  of  reference  12).  The  method- 
of-moment  solution  for  the  distributions  Iz(z) and  Jp{p)  is  described  in  reference  1.  These 
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Figure  1.  Monopole  Element  on  Disk  Ground  Plane  Above  Flat  Earth 
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currents  are  used  to  determine  the  far-zone  field  at  P(r,  8,0).  The  distribution  /Z(z)  and  Jp{p) 

include  the  contribution  of  the  surface  wave.  The  theory  which  follows  is  based  on  a  report  by 
Richmond  (see  Acknowledgments). 

Wherever  possible  the  notation  will  agree  with  that  of  references  1  and  10,  except  for  the 
following:  the  disk  radius  is  denoted  by  a  (instead  of  c);  the  monopole  element  radius  is  noted 
by  b  (instead  of  a);  the  disk  hole  radius,  equal  to  the  radius  of  coaxial  line  outer  conductor,  is 
denoted  by  b]  (instead  of  b);  the  Earth  complex  permittivity  is  denoted  by  e*  (instead  of  £2)*  die 
free-space  wave  impedance  is  denoted  by  T]  o  (instead  of  r\ );  and  the  origin  of  the  z  axis  is  on  the 
Earth's  surface  rather  than  on  the  disk  ground  plane.  These  changes  in  notation  are  made  to 
conform  with  the  notation  for  the  numerical  results  of  section  3.  The  time  dependence  is 
exp(jost).  The  parameters  of  free  space  are  denoted  by  (p0,e0),K=  co^p0£0,  and  rjG  =  ^p0/e0  . 

The  far-zone  electric  field  intensity  of  the  monopole/disk  antenna  may  be  regarded  as  the 
sum  of  the  field  E*  radiated  from  the  electric  currents  and  the  field  EM  radiated  from  the 
magnetic  frill  current  at  the  antenna  terminals.  To  calculate  these  fields  we  consider  the 
electric  current  density  J,  the  magnetic  current  density  M  radiating  in  free  space,  and  the  field 
reflected  from  the  air-Earth  interface.  In  these  calculations,  the  perfecdy  conducting  antenna 
structure  is  removed  and  replaced  with  the  equivalent  currents  J  and  M. 

One  successful  but  tedious  approach  to  the  far-zone  fields  starts  with  the  rigorous 
expressions  in  terms  of  Sommerfeld  integrals.  This  formulation  is  interpreted  as  a  plane- 
wave  expansion  that  includes  a  finite  spectrum  of  uniform  "space"  waves  plus  an  infinite 
spectrum  of  evanescent  plane  "surface"  waves.  Since  the  evanescent  waves  attenuate 
approximately  exponentially  with  their  height  above  the  Earth,  and  since  their  peak 
amplitude  relative  to  that  of  the  space  waves  approach  zero  with  increasing  distance  into  the 
far-field,  they  are  deleted  in  the  far-zone  field  derivations.  Finally,  the  method  of  stationary 
phase  is  applied  to  evaluate  the  remaining  integrals  asymptotically  as  the  observation  point 
recedes  to  infinity. 
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The  same  far- zone  field  expressions  can  be  derived  more  readily  via  Carson’s  reciprocity 
theorem  [22,23]  as  follows.  The  coordinate  origin  is  on  the  air-Earth  interface,  the  disk  in  the 
plane  z  =  z0,  and  the  wire  monopole  extends  from  z0  to  z0  +  h  on  the  z  axis.  A  point  in  the 
far-zone  region  is  specified  by  the  spherical  coordinates  (r,  6,  <p)  where  r  denotes  the  radial 
distance  from  the  origin  and  the  angle  9  is  measured  from  the  z  axis..  The  goal  is  to 
determine  the  far-zone  field  E(r,9,0)  in  the  free-space  region.  To  accomplish  this,  we  let 

denote  the  total  field  in  the  vicinity  of  the  origin  produced  by  an  infinitesimal 

electric  test  dipole  located  at  the  far-zone  point  (r,  0,  (p)  in  the  plane  0  =  0.  If  the  dipole  is 
oriented  in  the  6  direction  and  has  a  dipole  moment  p,  its  field  in  the  vicinity  (x,  y,  z)  of  the 
origin  is 


_ ,  _  -0jkr]opexp(-jkR' ) 

E0(x,y,z)  - - - - 

Anr 

R'  =  r  -  (x  sinO  +  z  cos 6) 

With  this  field  incident  on  the  air-Earth  interface,  the  total  magnetic  field  above  the  Earth  in 
the  vicinity  of  the  origin  is 

_ , .  yjkpexp{-jkr)  exp{jkxsin  0) 

H  (x,y,z)- - - - 

Anr 

[exp(  jkz  cos  9)  +  91  exp(-jkz  cos  6 )]  (3) 

where  51  denotes  the  plane- wave  Fresnel  reflection  coefficient  at  the  air-Earth  interface. 

The  Earth  conductivity,  permeability,  and  dielectric  constant  are  denoted  by  (<T,  pi,  £>). 
We  let  P2  =  A?-  in  which  case  the  reflection  coefficient  SR  (for  parallel  polarization)  is  given 
by 


(1) 

(2) 
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(4) 


9*  =  9l||  = 


(£*  l£0)COS  0  —  y(£*  I  £0)  —  Sin  6 
(£*  I £0) cos 6  +  J(£  *  l£a) -  sin 2  6 


The  complex  relative  permittivity  of  the  Earth  is 

(£*/£o)  =  £r-j<j/(0)£o)  (5 

From  equation  (3)  and  Maxwell's  curl  equations,  the  total  electric  field  intensity  above 
the  Earth  in  the  vicinity  of  the  origin  is 

-jkT]0pexp(-jkr)  exp(jkx  sin  0)  cos  6 

C,  =  —  — 

x  4nr 

[exp{jkz  cos  0)  -  *  exp(-  jkz  cos  0)]  (6 


Ely  =0 


£t  _  -jkrjoP  exp(-jkr)  exp(jkx  sind)  cos  6 
z  4nr 


[exp(jkz  cos  6)  +  *  exp(-  jkz  cos  0)] 


Carson’s  reciprocity  theorem  states  that 


\V  E dl  =  W^J  E'-M  H^ds 


where  E  is  the  total  far  zone  field  at  (r,  6,Q). 


On  the  right-hand  side  of  this  equation,  J  and  M  denote  surface-current  densities  on  the 
perfectly  conducting  monopole/disk  antenna  and  the  integration  extends  over  the  surface  of 
the  monopole/disk  antenna.  On  the  left  side  the  integration  extends  over  the  infinitesimal  test 
dipole  and  is  readily  evaluated  to  be  pEe(r,  8,  0),  so  the  reciprocity  theorem  reduces  to 

Ee{r, 8,4)  =  (1  /  p)  \\{J  ■  E‘  -  M  •  H‘ )  ds  (10) 

where  p  is  the  dipole  moment 


The  magnetic-frill  current  Af  is  given.  Upon  completion  of  the  moment-method  analysis, 

the  electric  current  density  J  is  known  on  the  vertical  wire  monopole  and  the  horizontal 

conducting  disk.  The  test  dipole  fields  are  given  above,  so  equation  (10)  contains  no 

unknown  quantities.  Thus,  evaluation  of  the  far- zone  field  Eq  (r,  8,  0)  of  the  monopole/disk 

antenna  is  now  simply  a  matter  of  performing  the  integrations  in  equation  (10).  If  we  start 
* 

with  a  0-oriented  test  dipole,  a  similar  analysis  shows  that  the  resulting  far-zone  0 
component  of  the  monopole/disk  antenna  is  Eq  (r,  8, 0)  =  0. 


2.2.  THE  FIELD  FROM  THE  MONOPOLE  ELEMENT 

In  far-zone  field  calculations  for  the  vertical  wire  monopole,  the  tubular  surface  current 
density  J  can  be  replaced  with  a  filamentary  line  source  I(z)  on  the  z  axis.  From  equations 
(8)  and  (10),  the  far-zone  field  from  the  vertical  wire  is  given  by 


Anr 

•P0  * I(z)[exp(jkzcos8)  +  yiecp(-jkz  cos8)]dz  (11) 
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The  vertical  wire  monopole  is  divided  into  L  segments  with  length  d'  =  L/h.  A  typical 
segment  (segment  l)  extends  from  zf  to  z|  on  the  z-axis,  with  the  following  current 
distribution: 


_  l{sink(zj-z)  +  ljsink(z-z{) 
Z  ~  sin  led' 


(12) 


The  current  entering  the  segment  at  the  bottom  is  l{  =  l{z[ ),  and  the  current  leaving  the 
segment  at  the  top  is  l{  =  /(z^). 

From  equations  (11)  and  (12),  the  far-zone  field  of  the  wire  monopole  is  given  by 

L 

E%(r,6,<f>)  =  C  ^ //  jerp(  jkz\  cos  6^-  A  exp{jkz[  cosOj 
/=i 

+9?  \exp{^-jkz[ cos d}-B exp^-jkz f  cos fljjj 
L 

+C  ^li{exp(jkz{  cos  8}-B  exp[  jkz{  cos  0) 

l 

+SR  \exp[^jkz[  cos  -  Aexp[- jkz^  cos  0)jj  (13) 

where 

A  =  cos  (kd)  +  j  cos  Osin  ( kd )  (14) 

B  =  cos(kd)  -  j  cos  0  sin  (kd)  (15) 
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c  jT)0  exp(-jkr ) 

Anr  sin(kd)  sin  6 


(16) 


On  the  lowest  wire  segment  (*  =  1) ,  the  current  at  the  bottom  is  /{  =IX.  On  the  highest 

segment  (t  =  L)  the  endpoint  currents  are  /{  =  1^  and  /£  =  0 ,  where  N  denotes  the  number 
of  equations  and  the  number  of  unknowns  in  the  moment-method  solution  for  the 
monopole/disk  antenna. 


23  THE  FIELD  FROM  THE  DISK  GROUND  PLANE 

The  electric  current  density  Jp{p)  on  the  perfectly  conducting  circular  disk  is  radially 
directed  and  independent  of  the  azimuthal  angle  0 .  The  disk  lies  in  the  plane  z  =  zo-  If 
(r  ,-0,>2o)  denotes  the  cylindrical  coordinates  of  a  source  point  on  the  disk,  the  far-zone  disk 
field  is  obtained  from  equations  (6),  (7),  and  (10)  as  follows: 

Eg  (i r ,  0,0)  =  ~jkVo  COS~-[exp(jkz0  cos  0)  -  91  exp(-jkz0  cos  0)]  (17) 

An  r 


\l\-nJP(p  ^cos(!>,'exP(jkP‘ cos  <P"sin  9)p'd<t>”dp' 

0"  =  0'-0  (18) 

Since  the  disk  current  density  Jp  is  independent  of  0',  one  integration  can  be  evaluated  as 
follows: 

f  cos 0  exp(jx  cos 0)  dp  =  2nj  J\(x)  (19) 

»  — IT 
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where  J\(x)  denotes  the  Bessel  function.  Beginning  at  this  point,  it  is  convenient  to  let  p 
(instead  of  p')  denote  the  radial  coordinate  of  a  source  point  on  the  disk.  From  equations 
(17)  and  (19),  the  far*zone  field  of  the  circular  disk  is 

=0.5  ki)0  \exp{- jkr)/r ]  cos  6 

[exp(jkz0 cos 0)  -  91  exp(-jkz0  cos  0)]  pJp(p)Jl(kpsin6)dp  (20) 

The  perfectly  conducting  circular  disk  is  divided  into  M  concentric  annual  zones.  A 
typical  zone  (zone  m)  has  an  inner  radius  pf\  an  outer  radius  pj1,  and  a  width 

b)/M .  Let  /”  denote  the  electric  current  entering  the  zone  at  pf,  and 
I 2  the  current  leaving  at  p-T .  Then  the  electric  surface  current  density  on  this  zone  is 


I\  sin  k(p2  -  p)  + 1 2  sin  k(p-p f* ) 
2  7i p  sin  kd 


(21) 


From  equations  (20)  and  (21),  the  far- zone  field  of  the  circular  disk  is  given  by 


E& (r, 6, <t>)  =  Vo ****(  6  [exp(jkz0 cos 0) -  91  exp{-jkz0 cos 0)] 

Anr  sin(kd)  1  J 

•  II  £'[/,"  *»  *(/>?  -  p)  +  >2  t(p  -  ft” )]  h  (*P  ««  0)  4¥)  (22) 

m=l  ^ 

On  the  first  zone  (m  =  1),  the  endpoint  currents  are  /{"  =  -Ix  and  /"  =  I 2-  0°  the  ^ast  zone 
(m  =  Af).  the  endpoint  currents  are  /r=/Af  and  Ij  =  0.  Numerical  integration  techniques 
are  required  in  evaluating  this  expression. 
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2.4  THE  FIELD  FROM  THE  MAGNETIC  FRILL 

The  perfectly  conducting  circular  disk  and  the  coaxial-fed  monopole  are  replaced  (via 
Schelkunoff  s  equivalence  principle)  with  equivalent  electric  and  magnetic  surface  currents 
radiating  in  free  space  over  the  flat  Earth.  The  equivalent  magnetic  surface-current  density, 
derived  in  section  2.4  of  reference  12,  is  given  by: 


-V/[p*n(fc,/b)],  b<p<b i 

0,  p  elsewhere 


(23) 


This  "magnetic  frill,”  located  at  z  =  z0.  is  centered  on  the  z-axis  and  has  inner  and  outer  radii 
of  b  and  bj,  respectively.  The  antenna  is  considered  to  be  transmitting,  with  a  voltage 
generator  (of  V  peak  volts)  at  the  terminals  and  the  coaxial  outer  conductor  at  zero  potential. 
The  free-space  field  of  the  magnetic  frill  is  analyzed  by  Tsai  [14,15].  From  equations  (3)  and 
(10),  the  far- zone  field  of  the  frill  is  given  by 

Eg  (r, 0, <p)  =  -r\exp{jkza cos 6)  +  91  exp(~jkza cos 0)1 

4  Jt  r 

•J;>  J  Jrij, (p')cos <p"exp( jkp' cos <p” sin 0)p'd<$>”dp'  (24) 

Since  the  magnetic  current  density  is  independent  of  <f>',  one  integration  can  be  performed 
with  the  aid  of  equations  (19)  and  (23)  to  obtain 


Ee  =  ~~ ^  - jgj- [exp(jkz0  cos  0)  +  91  exp(-jkz0  cos  0)]  •  J*  Jx  (kp  sin  0)  dp  (25) 


2-10 


The  final  integration  is  performed  as  follows: 


jjl(Px)dx  =  -J0px/P  (26) 

Thus,  the  field  of  the  magnetic  frill  is  given  by 


Ee  =  2rg/n(V/I)[eXp(^Z°  C°S  ^  +  ^  exp^~jkzo  cos  0)] 

•  [Ja(kb  sin  6)  -  J0(ka  sin  9)\/sin  6  (27) 

This  expression  can  be  simplified  with  the  following: 

J0(x)  =  l-x2/4,  x«l  (28) 

From  equations  (27)  and  (28),  the  far-zone  field  of  the  magnetic  frill  is  given  by 


Ee(r>  d,<t>) 


k2v[b2  -bf}exp(-jkr) 
8r  in(b  ~Jb) 


[eacp( jkzQ cos 6)  +  91  exp(-jkz0 cos 0)] sin 6,  kb x«\ 


(29) 


2.5  THE  TOTAL  FAR-ZONE  FIELD 

The  total  far-zone  field  E0(r,  6,<p)  defined  by  equation  (10)  in  the  free-space  region  is  the 
sum  of  the  fields  from  the  monopole  element,  the  disk  ground  plane,  and  the  magnetic  frill. 
Accordingly,  the  total  far-zone  field  in  the  free-space(air)  region  is  given  by 


2-11 


Ee(r,e,<p)=Ee(r,e)  =  E%(r,8)+Ei(r,d)  +  Elf(r,0) 


(30) 


where  Eq,Eq,  Eq  are  given  by  equations  (13),  (22),  and  (29),  respectively.  The  fields  £e, 
Eq,Eq,  and  Eq  are  uniform  with  azimuthal  angle  <p  because  of  the  azimuthal  symmetry  of 
the  antenna  geometry  in  figure  1. 


Consider  now  the  cases  where  the  Earth  medium  either  is  lossy  (<r  >  0)  or  is  free  space 
(or  =  0,  er  =  1).  The  total  far-zone  radiated  power  Pr  is  given  by 


j|£e(r,0)|2r2sm0  dd,  G>  0 


j  \Ee{r,6)f  r2  sinO  dd;  a- 0 


O 


(31) 


where  E0(r,0)  =  the  far-zone  field  (in  the  ffee-space  region)  given  by  equation  (30). 

T)0  =  (n0/e0)^2  =  free  space  wave  impedance  (ohms) 

For  the  case  of  a  >  0,  the  integrand  in  equation  (31)  is  integrated  over  only  the  hemisphere 
above  the  Earth  because  the  field  in  lossy  Earth,  relative  to  that  in  free  space,  approaches  zero 
at  large  radial  distances  r. 

The  antenna  directivity  d{6)  expressed  as  a  numeric  is  given  by 

d(6)=  2ffr2|£*M)|/(n*P,)  (32) 
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The  antenna  directivity  D(0),  expressed  in  decibels,  is  given  by 


D(0)  =  lOloglod(0)  ( dB )  (33) 

The  input  power  Pin  to  the  monopole  element  is  given  by 

^=(1/2)  Ke[V(0)/*(0)]  (34) 

where  V(0)  =  Peak  input  voltage  (volts).  The  input  voltage  V(0)  is  usually  set  equal  to 
1  volt  in  the  moment-method  analysis. 

1*(0)  =  Conjugate  of  the  peak  input  current  7(0)  at  the  base  of  the  monopole 
element.  This  current  is  solved  for  by  the  moment-method  analysis  in 
reference  1. 

The  input  impedance  Zin  is  given  by 

Zin~Rin  +  j  Xjn  =  V(0)//(0)  (35) 

where  Rin  and  Xm  are  the  input  resistance  and  reactance,  respectively. 

The  antenna  radiation  resistance  Rrad  is  defined  as 

K,*,  =  irr/\i(of  (36) 

The  antenna  radiation  efficiency  rj  is  defined  as 


(37) 
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For  the  case  of  firee-space  (a-0,er-  1 ),  the  radiation  efficiency  is  equal  to  unity  because 
the  monopole  element  and  the  disk  ground  plane  conductivities  are  assumed  to  be  infinite. 


SECTION  3 


NUMERICAL  RESULTS 


Numerical  evaluation  of  the  far- zone  field,  directivity,  radiation  resistance,  and  radiation 
efficiency  is  executed  by  Richmond's  computer  program  RICHMOND4  written  in 
FORTRAN  77,  with  double  precision  for  use  on  a  DEC  VAX  computer.  The  program 
RICHMOND4  uses  subroutines  from  Richmond's  computer  program  RICHMOND3  that 
determines  the  current  distributions  on  the  monopole  element  and  disk  ground  plane,  as  well 
as  the  input  current  I(Z0)  and  input  impedance  Z  =  V/I(z0).  Brief  descriptions,  listings,  and 
sample  outputs  by  Richmond  of  programs  RICHMOND3  and  RICHMOND4  are  given  in 
appendices  A  and  B,  respectively.  Programs  RICHMOND3  and  RICHMOND4  are 
extensions  of  programs  RICHMD1  and  RICHMD2,  respectively,  described  in  reference  12 
for  a  monopole  element  on  a  disk  ground  plane  in  free  space. 

Examples  of  numerical  results  are  presented  here  for  a  thin,  quarter-wave  monopole 
element  on  a  small  to  moderately  large  disk  ground  plane  resting  on  medium  dry  ground  at 
15  MHz  in  the  high-frequency  band  (b/X  =  10“*,  h/X  =  0.25,  Ina/X  =  0  to  8  wavenumbers, 
z0  -0,  £r-  15,  <r=  0.001  S/m,  tan5  =  60A oZer  =  0.08).  More  extensive  results,  in  the  form 
of  an  atlas  of  computer  plots,  are  presented  in  reference  21  as  a  function  of  Earth 
classification.  The  coaxial  line  feed  (bj/b  =  3.5)  has  a  negligible  effect  on  the  far-zone  field 
and  input  current  because  its  equivalent  magnetic  frill  of  outer  diameter  2 t\/X  (=  7  x  10*6 
wavelengths)  has  a  radiation  resistance  that  is  small  compared  to  that  of  the  monopole 
element  of  length  h/X  (=  0.25  wavelengths).  In  the  numerical  results,  the  monopole  element 
was  divided  into  four  segments.  The  disk  was  segmented  into  equal- width  annular  zones, 
whose  numbers  varied  from  seven  for  ka  -  0.025, 0.25, 0.50;  16  for  ka  =  0.75  through  5.25; 
17  for  ka  =  5.5;  18  for  ka  =  5.75  and  6.00;  19  for  ka  =  6.25;  20  for  ka  =  6.5;  21  for  ka  =  6.75, 
7.0;  22  for  ka  -  7.25;  23  for  ka  -  7.50;  and  24  for  ka  =  7.75  and  8.0.  Results  are  compared 
with  those  for  a  perfect  ground  plane  (£r  =  1.0,  a  =  °°)  and  for  an  Earth  permittivity  equal  to 
that  of  free  space  (£r  =  1.0,  c  =  0).  The  results  for  perfect  ground,  medium  dry  ground,  and 
free  space  are  identified  in  the  following  figures  as  Case  1,  Case  5,  and  Case  11,  respectively. 
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The  elevation  numeric  directivity  patterns  for  disk  radii  2  tw/A  =  0.025,  3.0, 4.0,  5.0,  and 
6.5  wavenumbers  are  shown  as  polar  plots  on  the  same  linear  scale  in  figures  2  through  6, 
respectively.  In  the  presence  of  Earth  (Case  5),  the  directivity  patterns  are  approximately 
independent  of  disk  radius.  The  Earth  softens  the  edge  of  the  ground  plane  and  minimizes 
changes  in  directive  gain  resulting  from  ground  plane  edge  diffraction.  The  peak  directivity 
(see  figure  7)  is  within  0.5  dBi  of  that  for  a  perfect  ground  plane.  The  direction  of  peak 
directivity  (see  figure  8)  is  approximately  30°  above  the  horizon  with  variations  of  less  than 
4°  for  ground  plane  radii  0  <  Inajk  <  8  wavenumbers.  The  directivity  at  angles  of  incidence 
near  the  horizon  (see  figures  9  through  13)  for  0  Ina/X  <,  8  wavenumbers  has  no 
improvement  over  that  with  no  ground  plane  at  all  and,  in  fact,  decreases  periodically  with 
increasing  disk  radius  by  as  much  as  1  dB.  The  directivity  at  angles  of  incidence  of  82°,  84°, 
86°,  88°,  and  90°  are  approximately  4  dB,  5  dB,  7  dB,  13  dB,  and  «» dB,  respectively,  below 
the  peak  directivity  for  these  disk  radii. 

The  directivity  on  the  horizon  (see  figure  13)  is  -«>  dB  because  of  the  space  wave 
multipath  null  for  Earth  surface  reflection  at  a  grazing  angle  of  0°.  In  actuality,  the  field  on 
the  radio  horizon  is  not  zero  because  of  the  leaky  evanescent  surface  wave  that  is  generated 
in  the  air  medium  in  proximity  to  the  air-Earth  interface  [13].  The  surface  wave  has  an 
evanescent  field  in  the  air-medium  only,  but  leaks  energy  into  the  Earth  medium,  not  into  the 
air  medium.  The  amplitude  of  the  space  wave  in  the  direction  of  peak  directivity  approaches 
zero  with  increasing  distance  into  the  far-zone. 

The  theoretical  numeric  directive  gain  of  electrically  short  monopole  elements  on  ground 
planes  resting  on  lossy  Earth  may  be  approximated  by  an  expression  of  the  form  [13] 


dr{0)  = 


\Acosm  0sinn  6;  O£0£tc/2  rad,  m>0,n>\ 
[0,  —  zr/2  £  6  <  0  rad 


(37) 
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f  =  15  MHz 

h/k=0.25,  b/X.=1.0  x  10  6,  Vk=° 

Case  1,  Perfect  Ground  (£r=1.0,  o-°°) 
Case  5,  Medium  Dry  Ground 
(e  =15.0,  o=0. 001  S/m) 

Case  11,  Free  Space  (e=1.0,a=0) 


/  \ 


case  5  / 


?N-,  \  ,-A 

vOA.  case  11.--''  \  I 


16S*-  ISO*  135*  V 


Figure  6.  Numeric  Directive  Gain  Polar  Plot,  2  jeo/A  =  6.5 
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case  1 
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Figure  II.  Directive  Gain  at  Four  Degrees  above  the  Horizon 
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Figure  12.  Directive  Gain  at  Two  Degrees  above  the  Horizon 


The  exponents  m  and  n  are  chosen  to  yield  a  peak  directivity  in  a  desired  direction,  and  a  null 
at  0  =  0  and  re/2  radians.  The  coefficient  A  is  chosen  to  satisfy  the  condition 


In  m2 

(l/4/r)J  jdr(0)sin0  d<f>  =  1 
o  o 

Accordingly, 


A 


Odd 


The  directivity  of  equation  (37)  has  a  null  in  the  direction  of  zenith,  and  the  horizon  has  a 
peak  directivity  comparable  to  that  for  a  perfect  ground  plane. 

An  analytical  expression  that  approximates  the  directivity  obtained  by  numerical  methods 
for  medium  dry  ground  and  2tcaA  =  3  (see  figure  14)  is  given  by 


dr(0)  = 


\\Ocos0sir?  0;  O£0<n/2  rad 
[0,  -tt/2 £0<O rad 


(38) 


In  the  absence  of  Earth  (Case  1 1),  the  directivity  patterns  (see  figures  2  through  6)  are 
strong  functions  of  the  disk  radius  because  ground-plane  edge  diffraction  is  more 
pronounced.  The  peak  directivity  (see  figure  7)  varies  from  approximately  2  dBi  to  5  dBi. 
The  angle  of  incidence  of  peak  directivity  (see  figure  8)  varies  from  0°  to  32°.  The  large 
changes  in  angle  of  peak  directivity  at  2 iuil X  =  5.5  wavenumbers  do  not  represent  significant 
changes  in  peak  directivity  because  of  the  broad  3-dB  beam  width  of  the  directivity  pattern. 
The  jump  in  angle  of  peak  directivity  between  2na/X  =  5.5  and  5.75  wavenumbers 
corresponds  to  a  change  in  beamshape  (compare  figures  5  and  6).  The  directivity  on  the 
horizon  (see  figure  13)  varies  from  1.88  dBi  for  2na!\  =  0  wavenumbers  to  the  asymptotic 
value  of  -0.88  dBi  for  large  disk  ground  planes  of  finite  radius. 
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Numeric  Directive  Gain,  df  (0) 


ha  =  025,  b/X  =  1 .0  x  1 0'6,  z<A  =  0, 1  =  15  MHz 
Case  5,  Medium  Dry  Ground 
(£r  =  15.0,  a  =  0.001  S/m) 

Figure  14.  Numeric  Directive  Gain  of  a  Quarter-Wave  Element  on  a  Disk  Ground 
Plane  Resting  on  Medium  Dry  Ground,  2xa/X  =  3.0 
(a)  Richmond's  Method-of-Moments;  (b)  10  cos  0  sin3  6 
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The  radiation  resistance  (see  figure  15)  increases  aperiodically  with  increasing  disk 
radius.  The  aperiodicity  is  more  apparent  in  the  absence  of  Earth  because  ground-plane  edge 
diffraction  is  more  pronounced.  The  radiation  efficiency  (see  figure  16)  increases 
monotonically  with  increasing  disk  radius  in  the  presence  of  Earth:  from  0.21  for  2rtaA  =  0, 
to  0.69  for  2itaA  =  8,  and  to  1.0  for  2;taA  =  In  the  absence  of  Earth,  the  radiation 
efficiency  is  equal  to  unity  because  the  monopole  element  and  disk  are  assumed  to  be  of 
infinite  conductivity.  The  reason  why  the  radiation  efficiency  is  so  small  for  small  disk 
ground  planes  in  close  proximity  to  Earth,  regardless  of  whether  the  Earth  is  lossy  (o  >  0)  or 
is  a  pure  dielectric  (o=  0),  is  because  most  of  the  available  input  energy  is  directed  into  the 
Earth  by  the  leaky  evanescent  surface  wave  generated  by  the  spherical  wave  source  (the 
monopole  element)  in  the  air  medium  in  proximity  to  the  air-Earth  interface  [13]. 

Richmond's  method -of-moments  model,  in  solving  for  the  input  current  l(zQ)  indirectly, 

includes  the  Si  rface  wave  and  its  affect  on  the  far- zone  radiation  resistance  and  radiation 
efficiency.  Although  this  paper  is  restricted  to  the  calculation  of  the  far-zone  field  above  the 
Earth,  Richmond's  moment  method  analysis  for  the  element  and  disk  current  distributions  can 
also  be  used  to  calculate  the  near-zone  field  including  that  of  the  surface  wave.  This  latter 
effort  has  not  yet  been  undertaken. 
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SECTION  4 

VALIDATION  OF  NUMERICAL  RESULTS 


Several  approaches  have  been  used  in  validating  the  numerical  results  from  the  programs 
RICHMOND 3  and  RICHMOND4.  These  approaches  include  validation  by  comparison  with 
results  from  the  limiting  case  of  disk  ground  planes  in  free  space;  the  limiting  case  of  ground 
planes  of  zero  extent  in  proximity  to  Earth;  the  limiting  case  of  a  perfect  ground  plane  of 
infinite  extent;  Wait-Surtees  model  for  input  impedance;  Wait-Waiters  model  for  gain;  and 
the  Numerical  Electromagnetics  Code  (NEC)  for  radiation  efficiency. 

4.1  LIMITING  CASE  OF  DISK  GROUND  PLANES  IN  FREE  SPACE 

In  the  limiting  case  of  disk  ground  planes  in  free  space,  numerical  results  from  programs 
RICHMOND3  and  RICHMOND4  agree  with  results  from  programs  RICHMD1 
and  RICHMD2.  The  method-of-moments  programs  RICHMD1  and  RICHMD2,  for  a 
monopole  element  on  a  disk  ground  plane  in  free  space,  have  received  extensive  validation 
[12].  In  reference  12,  numerical  results  for  electrically  thin  monopole  elements  were 
compared  with  results  from  Brillouin-Stratton  induced  electromotive  force  (EMF)  method  for 
ground  planes  of  zero  extent;  Bardeen's  integral  equation  method  for  ground-plane  radii 
0<ka<  2.75  wavenumbers;  Leitner-Spcnce  method  of  oblate  spheroidal  wave  functions  for 
ground  plane  radii  3.0  £  ka  <  6.5  wavenumbers;  Awadalla-McClean  moment  method 
combined  with  the  geometric  theory  of  diffraction  for  ground-plane  radii  8.5  <  ka  <  <» 
wavenumbers;  and  the  method  of  images  for  ka  =  <».  Consistent  and  excellent  agreements  of 
results  were  achieved  by  the  RICHMD1  and  RICHMD  2  programs. 


4.2  LIMITING  CASE  OF  GROUND  PLANES  OF  ZERO  EXTENT 

In  the  limiting  case  of  ground  planes  of  zero  extent  in  proximity  to  Earth,  program 
RICHMOND4  results  for  the  directivity  of  a  quarter-wave  monopole  element  with  a  disk 
ground  plane  of  radius  ka  =  0.025  wavenumber  resting  on  medium  dry  Earth  (see  Case  5  of 
figure  2)  were  compared  with  results  for  ka  =  0  from  a  Fresnel  reflection  model  (MITRE 
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Program  MODIFIED  IMAGES)  and  Lawrence  Livermore  Laboratory’s  method-of-moments 
program  NEC-3  using  the  Sommerfeld  option.  Programs  RICHMOND4,  MODIFIED 
IMAGES,  and  NEC-3  gave  identical  directivity  patterns  with  absolute  values  of  directivity 
that  agreed  to  within  0.04  dBi.  The  reason  for  the  close  agreement  is  that  the  directivity 
does  not  depend  upon  the  absolute  accuracy  of  the  antenna  input  current. 

Radiation  resistance  and  radiation  efficiency  do  depend  upon  the  absolute  accuracy  of  the 
antenna  input  current.  RICHMOND4  results  of  radiation  resistance  and  radiation  efficiency, 
for  the  above  case  and  various  types  of  Earth,  are  compared  in  table  1  with  results  from 
NEC-3  (but  not  MODIFIED  IMAGES  because  the  omission  of  the  surface  wave  in  the 
Fresnel  coefficient  model  affects  the  radiation  efficiency  and  radiation  resistance,  but  not 
directivity).  The  results  differ  by  approximately  10%  for  radiation  resistance  and  by  more 
than  25%  for  radiation  efficiency.  These  differences  are  attributable  to  the  difference  in 
charge  density  at  the  base  of  the  monopole  element  by  a  factor  of  4000  resulting  from  the 
different  configurations  of  the  two  models  [16].  In  NEC-3,  the  current  produced  by  the 
charge  distribution  is  discharged  into  the  Earth  through  an  element  of  radius  10'6 
wavelengths,  whereas  in  RICHMOND4  the  current  is  discharged  into  the  Earth  through  a 
ground  plane  of  radius  4  x  10' 3  wavelengths.  The  NEC-3  results  for  the  radiation  efficiency 
of  a  quarter- wave  monopole  element  is  augmented  by  a  128-radial-wire  ground  plane  of 
radius  0.01  wavelengths  (see  section  4.6).  An  increase  in  the  number  of  monopole  segments 
from  4  to  20  in  RICHMOND4  has  no  significant  effect  in  modifying  the  table  7  results  for 
radiation  efficiency. 

4.3  LIMITING  CASE  OF  A  GROUND  PLANE  OF  INFINITE  EXTENT 

In  the  limiting  case  of  a  perfect  ground  plane  of  infinite  extent,  the  monopole  element  of 
length  h  may  be  modeled  by  the  method-of-images  as  a  free-space  dipole  of  half-length  h,  but 
with  twice  the  dipole  input  current,  one-half  the  dipole  impedance,  twice  the  dipole 
directivity  in  the  upper  hemisphere,  and  zero  times  the  dipole  directivity  in  the  lower 
hemisphere.  Richmond  has  written  a  program,  RICHMD6,  that  uses  a  sinusoidal-Galerkin 
method  of  moments  to  compute  the  input  impedance,  current  distribution,  and  far- zone  field 
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Table  1.  Radiation  Resistance  and  Efficiency  of  a  Vertical,  Quarter-Wave, 
Monopole  Element  on  Flat  Earth;  f  =  15  MHz,  b fk  =  1.0  x  10"^ 


Earth 

Classification 
(er,  a  S/m) 

|  Radiation  Resistance  (Ohms) 

|  Radiation  Efficiency  (Numeric) 

•NEC-3 

♦♦RICHMD4 

ncffljijs 

ImTmniT?! 

•NEC-3 

••RICHMD4 

•••Percent 

Difference 

Seawater 

(70,5) 

34.0 

29.5 

15.0 

0.823 

0.799 

29.4 

Fresh  water 
(80,  3.0  x  10-2) 

19.1 

17.3 

10.4 

0.273 

0.347 

34.3 

Wet  ground 
(30,  1.0  x  lO-2) 

14.5 

13.2 

10.3 

0.144 

0.229 

36.9 

Medium  dry 
ground 

(15,  1.0  x  10"3) 

11.5 

10.5 

10.3 

0.163 

0.210 

22.2 

Very  dry  ground 
(3,  1.0  x  1(H) 

6.2 

5.7 

9.6 

0.091 

0.145 

37.6 

19.1 

17.3 

9.4 

0.375 

0.378 

0.8 

Ice(-1°Q 
(3.  9.0  x  10-5) 

6.2 

5.7 

9.6 

0.091 

0.148 

38.8 

Ice  (-10°Q 
(3,  2.7  x  lO"5) 

6.2 

5.7 

9.5 

0.136 

0.171 

20.8 

Average  land 
(10,  5.0  x  10"3) 

9.9 

9.0 

10.3 

0.044 

0.105 

58.3 

*  Number  of  element  segments,  N  =  25;  voltage  source  excitation  at  N  =  1 
**  Disk  ground  plane  radius,  2naA  =  0.025  wavenumbers 
***  KNEC-3  -  RICHMD4)/  RICHMD4I  x  100 


4-3 


of  the  equivalent  free-space  dipole.  A  listing  of  program  RICHMD6  is  given  in  appendix  C. 
Numerical  results  for  input  impedance  are  in  reasonable  agreement  with  those  from  King- 
Middleton  theory  [17].  For  example,  for  htX  =  0.25  (corresponding  to  kh  =  it  12)  and 
hlb  -  16.56  (corresponding  to  ft  =  7),  RICHMD6  results  for  the  monopole  input  impedance 
are  Zin  =  46.52  +  j  15.97  ohms  which  differ  from  the  King-Middleton  results  of  Z,n  =  47.85  + 
j  18.50  ohms  by  2.8  and  13.7%  for  input  resistance  and  input  reactance,  respectively. 
RICHMD6  results  for  directivity  are  almost  identical  to  the  well-known  results  for  a  thin, 
quarter- wave  monopole  on  a  perfect  ground  plane  [12]. 


4.4  COMPARISON  WITH  WAIT-SURTEES  MODEL  FOR  INPUT  IMPEDANCE 

Program  RICHMOND4  results  for  the  input  impedance  of  a  monopole  element  with  a 
disk  ground  plane  resting  on  flat  Earth  have  been  compared  by  Richmond  [1]  with  those 
obtained  from  a  Wait-Surtees  model  [18].  In  reference  1,  the  Wait-Surtees  results  for  input 
reactance  are  inadvertently  given  for  a  disk  ground  plane  in  free  space  rather  than  for  a  disk 
ground  plane  on  flat  Earth.  RICHMOND4  results  for  input  resistance  and  input  reactance 
are  compared  in  figures  17  and  18,  respectively,  with  those  obtained  from  a  program  WAIT- 
SURTEES  written  by  Richmond  and  based  on  the  Wait-Surtees  model.  Program  WAIT- 
SURTEES,  described  in  Appendix  D,  incorporates  results  from  program  RICHMD6  for  the 
input  impedance  of  a  monopole  element  on  a  perfect  ground  plane.  The  RICHMOND4 
results  are  in  close  agreement  with  WAIT-SURTEES  results,  except  at  small  ground-plane 
radii  less  than  approximately  to  =  1.0  wavenumber  for  which  the  Wait-Surtees  model  is  not 
accurate.  Richmond  [10]  has  compared  RICHMD1  results  with  WAIT-SURTEES  results  for 
the  input  impedance  of  a  monopole  element  on  a  disk  ground  plane  in  free  space  and 
obtained  similar  agreement  as  above,  but  for  ground-plane  radii  greater  than  approximately 
to  =  2.0  wavenumbers.  The  RICHMOND4  results  in  figure  18  for  input  reactance  should  not 
have  a  local  minimum  at  to=  0.75.  A  nonconvergent  result  was  obtained  at  to  =  0.75 
because  of  over-segmentation  of  the  disk  when  the  number  of  disk  annular  zones  was 
abruptly  increased  from  seven  at  to-  0.5  to  sixteen  at  to-  0.75. 
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Figure  17.  Input  Resistance 
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Figure  18.  Input  Reactance 


4.5  COMPARISON  WITH  WAIT- WALTERS  MODEL  FOR  GAIN 


Numerical  results  of  directivity  and  radiation  efficiency  from  Richmond's  method-of- 
moments  program  RICHMOND4  cannot  be  validated  against  models  based  on  Monteath's 
compensation  theorem  [5-8,19]  or  Sommerfeld's  attenuation  function  [9]  because  those 
models  yield  only  the  gain  (the  product  of  directivity  and  radiation  efficiency)  rather  than 
directivity  and  radiation  efficiency  as  separate  entities.  Nevertheless,  it  is  of  interest  to 
compare  RICHMOND4  results  for  gain  with  those  from  the  Wait-Waiters  model  [6,7,8,19] 
based  on  Monteath's  compensation  theorem. 

First  consider  the  Wait-Waiters  model.  The  gain  G  (ka,  if/)  -  G  (0,  if/)  (dB)  of  an 
electrically  short  monopole  element  on  a  disk  ground  plane  with  radius  ka  wavenumbers 
relative  to  that  without  a  disk  ground  plane  ( ka  =  0)  is  shown  in  figure  2  of  reference  8  and  in 
figure  23.26  of  reference  19  for  ka  =  10,  er  =  9,  and  <7=0.  The  Wait-Waiters  model  of 
reference  8  computes  the  magnetic  field  intensity  H  (ka,  if/ )  with  a  disk  ground  plane  as  a 
function  of  the  grazing  angle  if/  (the  complement  of  the  angle  of  incidence  0)  relative  to  that 
with  no  ground  plane.  At  a  grazing  angle  if/  =  2°,  the  Wait-Waiters  model  gives  a  relative 
gain  of  G(10, 2)  -  G  (0, 2)  =  4.5  dB 

Now  consider  the  Richmond  model.  Program  RICHMOND4  results  for  a  quarter-wave 
monopole  element  on  a  disk  ground  plane  of  radius  ka  =  8  wavenumbers  on  medium  dry 
ground  (er  =  15.0,  <7  =  0.001  S/m)  gives  a  directivity  at  a  grazing  angle  if/  =  2°,  of  D(8, 2) 

=  -8.6  dBi  (see  figure  12)  and  a  radiation  efficiency  tj  =  0.69  =  -1.6  dB  (see  figure  16).  The 
gain  G  (8, 2)  =  -8.6  dBi  -1.6  dB  =  -10.2  dB;  for  ka  =  0  and  if/  =  2°,  D( 0,2)  =  -7.9  dBi  (see 
figure  12)  and  the  radiation  efficiency  Tf  =  0.21  =  -6.8  dB  (see  figure  16).  The  gain  G(0,2)  = 
-7.9  dBi  -  6.8  dBi  =  -14.7  dBi.  The  relative  gain  G  (8,2)  -  G(0,2)  =  -10.2  dBi  +  14.7  dBi  = 
4.5  dB. 

The  RICHMOND4  and  Wait-Waiters  results  of  4.5  dB  for  relative  gain  are  identical  for 
these  similar  cases. 
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4.6  COMPARISON  WITH  NEC  FOR  RADIATION  EFFICIENCY 

Numerical  results  of  radiation  efficiency  obtained  from  programs  RICHMOND4,  NEC-3, 
and  NEC-GS  are  compared  in  figure  19  for  the  radiation  efficiency  of  a  quarter- wave 
monopole  element  with  small  ground  planes  on  or  just  above  medium  dry  Earth  as  a  function 
of  the  ground-plane  radius.  RICHMOND4  results  are  for  disk  ground  planes  (see  figure  16). 
NEC-3  results  are  for  a  ground  plane  of  zero  extent  (see  table  1).  NEC-GS  results  are  for 

radial- wire  ground  planes  whose  wires  have  a  radius  bw  =  10-5  wavelengths  [16,20],  The 
results  for  disk  ground  planes  are  in  close  agreement  with  those  for  ground  planes  with  128 
radial  wires. 
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Normalized  groundscreen  radius,  2na/X  (wave  number) 

Figure  19.  Radiation  Efficiency  of  a  Quarter-Wave  Monopole  Element  with 
Different  (Zero-Extent,  Radial-Wire,  and  Disk)  Ground  Planes 
on  or  just  above  Medium  Dry  Earth 


SECTION  5 
CONCLUSIONS 


Richmond's  moment-method  results,  for  the  current  distribution  and  input  impedance  of  a 
monopole  element  on  a  disk  ground  plane  above  flat  Earth,  are  used  to  obtain  the  far-zone 
field,  directivity  pattern,  radiation  resistance,  and  radiation  efficiency.  This  model  for  a  disk 
ground  plane  complements  the  NEC  method-of-moments  model  of  Burke,  et  al.  for  a  radial- 
wire  ground  plane. 

Method-of-moments  models,  unlike  models  based  on  Sommerfeld's  attenuation  function 
or  variational  models  based  on  Monteath's  compensation  theorem,  determine  the  directivity 
and  radiation  efficiency  as  separate  entities  rather  than  lumping  them  together  as  a  product  to 
yield  the  antenna  gain.  Other  advantages  of  the  method-of-moments  models  are  more  exact 
determination  of  current  distributions;  applicability  to  electrically  small  ground  planes;  direct 
determination  of  ground-plane  edge  diffraction;  and  avoidance  of  analytical  restrictions  on 
evaluating  Sommerfeld's  integral.  The  segmentation  of  ground  planes  in  method-of-moments 
models  restricts  the  models  to  ground  planes  that  are  sufficiently  small  so  that  computer 
computational  capacity  and  precision  are  not  exceeded. 

The  far-zone  field  in  the  free-space  (air)  region  is  determined  as  the  sum  of  direct  and 
indirect  (reflected  from  the  Earth)  fields  from  the  monopole  element,  disk  ground  plane,  and 
the  magnetic  frill  of  the  coaxial-line  feed  excitation.  The  far-zone  direct  fields  from  the 
monopole  element  and  disk  ground  plane  are  determined  from  the  method-of-moments 
solution  for  their  current  distributions.  The  far-zone  indirect  fields  are  determined  using  the 
plane-wave  Fresnel  reflection  coefficient.  The  significant  contribution  of  the  surface  wave  to 
the  far-zone  field  at  or  near  the  air-Earth  interface  is  not  considered,  but  is  small  compared  to 
the  far-zone  field  in  the  direction  of  peak  directivity.  However,  the  significant  effect  of  the 
surface-wave  in  determining  input  current  and  radiation  efficiency  are  included  in  the  present 
analysis. 
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Examples  of  numerical  results  are  presented  for  the  directivity  pattern,  peak  directivity, 
radiation  resistance,  and  radiation  efficiency  of  a  thin,  quarter- wave  monopclc  clement  on 
small  to  moderately  large  disk  ground  planes  (of  radius  0  to  8  wavenumbers)  resting  on 
medium  dry  Earth.  Results  are  compared  with  these  for  a  ground  plane  of  infmite  extent  and 
for  ground  planes  in  free  space.  In  the  presence  of  Earth,  the  directivny  patterns  are 
approximately  independent  of  disk  radius  for  ground-plane  radii  at  least  as  large  as  eight 
wavenumbers.  The  peak  directivity  is  within  0.5  dBi  of  that  for  a  perfect  ground  plane.  The 
direction  of  peak  directivity  is  approximately  30°  above  the  horizon.  The  directivity  at 
angles  of  incidence  of  82°,  84°,  86°,  88°,  and  90°  are  approximately  4  dB,  5  dB,  7  dB,  13  dB, 
and  «>  dB,  respectively,  below  the  peak  directivity.  The  numeric  directivity  is  given 
approximately  by  the  empirical  expression  10  cos  0  sin^  0  in  the  hemisphere  above  the  Earth 
and  by  zero  in  the  hemisphere  below  the  Earth.  The  radiation  efficiency  increases 
monotonically  with  increasing  disk  radius  in  the  presence  of  Earth:  from  0.21  for  a  ground 
plane  of  zero  extent  to  0.69  for  a  ground-plane  radius  of  eight  wavenumbers. 

Numerical  results  from  Richmond’s  method-of-moments  computer  programs 
RICHMOND3  and  RICHMOND4  for  a  monopole  element  with  a  disk  ground  plane  above 
flat  Earth  are  in  good  agreement  with  results  known  from  other  models  in  the  limiting  cases 
of  disk  ground  planes  in  free  space,  disk  ground  planes  of  zero  extent  in  proximity  to  Earth, 
and  a  perfect  ground  plane.  RICHMOND3  results  for  input  impedance  are  in  good 
agreement  with  results  from  a  Wait-Surtees  variational  model,  except  for  ground-plane  radii 
less  than  approximately  one  wavenumber  for  which  the  Wait-Surtees  model  is  not  accurate. 
A  RICHMOND4  result  for  antenna  gain  is  in  agreement  with  a  result  from  a  Wait-Waiters 
variational  model.  RICHMOND4  results  of  radiation  efficiency  are  in  close  agreement  with 
NEC-GS  method-of-moment  results  for  ground  planes  with  a  large  number  of  radial  wires. 
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COMPUTER  PROGRAM  RICHMOND3  FOR  THE  INPUT  IMPEDANCE  AND 
CURRENT  DISTRIBUTION  OF  A  MONOPOLE  ELEMENT  ON  A  DISK  GROUND 

PLANE  ABOVE  FLAT  EARTH 


COMPUTER  PROGRAM  RICHMOND3 
MONOPOLE  ANTENNA  ON  CIRCULAR  DISK  OVER  FLAT 

EARTH 

(Current  Distribution  and  Impedance) 
by 

Jack  H.  Richmond 
January  29,  1990 


INTRODUCTION1 

Appendix  I  presents  the  computer  program  RICHMOND3  together  with 
all  the  necessary  subroutines.  This  FORTRAN  program  calculates  the 
current  distribution  and  impedance  of  a  monopole  antenna  mounted  at  the 
center  of  a  circular  disk  over  the  flat  lossy  earth. 

See:  [J.  H.  Richmond,  “Monopole  Antenna  on  Circular  Disk  Over  Flat 
Earth,”  IEEE  Transactions,  Vol.  AP-33,  pp.  633-637,  June  1985.] 

Comment  statements  have  been  inserted  in  the  main  program  and  the 
subroutines  to  assist  the  user.  Only  a  few  brief  comments  will  be  required 
in  this  Introduction. 

RICHMOND3  performs  all  calculations  with  double  precision.  The  the¬ 
oretical  basis  for  this  program  is  presented  in  the  above  published  paper, 
and  the  notation  in  the  program  corresponds  with  the  notation  in  the  pa¬ 
per  with  one  exception.  The  outer  radius  of  the  disk  is  denoted  by  b  in  the 
program  and  by  c  in  the  paper. 

In  Appendix  I,  Table  I  presents  the  antenna  impedance  as  calculated 
with  RICHMOND3  on  a  VAX  computer  for  the  following  disk  radii:  BL 
=  0.1,  0.2,  0.3  and  0.4.  The  antenna  impedance  (in  free  space  and  on  a 
lossy  flat  earth)  agree  closely  with  the  original  calculations  obtained  on 
a  DATA  CRAFT  computer  in  May  1979.  Table  II  presents  the  current 
distributions  on  the  monopole  and  the  disk  (in  free  space  and  on  a  lossy 

1  Appreciation  is  expressed  to  The  MITRE  Corporation  for  sponsoring  this  report . 

The  computer  program  RICHMOND3  was  developed  (in  single  precision)  in  1979  with 
other  sponsorship. 
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flat  earth)  with  BL  =  0.1.  These  results  also  agree  closely  with  the  original 
calculations  of  May  1979. 

Table  III  presents  the  antenna  impedance  (in  free  space  and  on  a  fiat 
earth)  as  calculated  with  a  VAX  with  single  precision.  Comparison  with 
Table  I  indicates  that  the  need  for  double  precision  is  marginal  for  this  case. 

In  the  original  program,  subroutine  CROUT  was  employed  to  solve 
the  simultaneous  linear  equations.  In  RICHMOND3,  CROUT  U  replaced 
with  CMINV  which  employs  full  pivoting  (on  rows  and  columns)  whereas 
CROUT  does  not  pivot.  On  the  other  hand  CMINV  i6  presumably  slower 
(in  solving  large  matrix  equations)  because  it  inverts  the  matrix,  whereas 
CROUT  solves  the  equations  without  inverting. 

The  current  distribution  will  be  printed  from  CMINV  if  IWCJ  =  1,  but 
the  printout  will  be  suppressed  if  IWCJ  =  0. 

Diagnostic  data  will  be  printed  from  several  subroutines  if  IWZ  =  1. 
This  printout  i6  suppressed  if  IWZ  =  0. 

The  integer  NPH  controls  the  numerical  integrations  in  several  subrou¬ 
tines.  NPH  determines  the  number  of  times  the  integrand  is  to  be  sampled 
with  Simpson's  rule.  The  value  NPH  =  6  usually  gives  a  suitable  com¬ 
promise  between  accuracy  and  computational  expense.  A  larger  value  will 
increase  the  expense,  and  it  may  improve  the  accuracy  in  some  cases. 

TL  denotes  the  thickness  of  the  circular  disk,  measured  in  free-space 
wavelengths.  To  promote  convergence  of  the  moment  method  (as  NEQ  i6 
increased),  the  value  TL  =  AL/100  is  recommended  regardless  of  the  true 
thickness  of  the  metallic  circular  disk.  This  result  is  rather  unexpected,  and 
the  interpretation  is  not  totally  understood.  Of  course,  it  is  assumed  that 
the  disk  thickness  TL  (as  well  as  the  monopole  wire  radius  AL)  is  much 
smaller  than  the  wavelength. 
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Appendix  I.  R3EHHOM>3  Program  Listing 


C  RICHMOND? 

C  MONOPOLE  AT  CENTER  or  CIRCULAR  DISK  OVER  FLAT  EARTH.  R1OM0ND3. 1 

C  DOUBLE  PRECISION. 

C  CURRENT  DISTRIBUTION  AND  IMPEDANCE. 

C  SEE:  RICHMOND,  "MONOPOLE  ANTENNA  ON  CIRCULAR  DISK  OVER  FLAT  EARTH", 

C  IEEE  TRANS.,  VOL.  AP-33,  PP .  633-637,  JUNE  1985. 

C  LINK:  BES10,CISI,CMINV,DZ11,DZDD,DZWD,DZWW,EXPJ,  GRILL, 

C  ODD,  QDM,  OHM,  SKEW,  SKEWS,  SKEHT,  SPART,  ZSDM,  ZSMM 

IMPLICIT  REAL* 8  (A-H) ,  (P-Z) 

COMPLEX* 16  CJ (30) ,VJ (30) , ZJ (30) ,VIJ(30, 30) , ZIJ (30, 30) 

COMPLEX* 16  Til , DET, EC, D1 1 , D12 , D21 , D22 , DZ1 J, DV1 , Ml 1 
COMPLEX* 16  P11,P12,P21,P22, ZDD, ZDM, ZMD, ZMM, Zll, ZD, Z22, Z12, Z21 
DIMENSION  FB  (500)  ,LLL  (30)  ,  MMM  (30) 

DATA  EO,UO/8.85418S33677E-12, 1 . 25663706144E-6/ 

DATA  ETA, PI, TP/376. 730366239,  3. 14159265359, 6.28318530718/ 

DATA  ICC, IFB/30,  500/ 

1  FORMAT (IX, 215, 7E15. 4) 

2  FORMAT  (IX,  7F17 . 8) 

S  FORMAT  (1H0) 

C  AL  -  RADIUS  OF  MIRE  IN  WAVELENGTHS. 

C  HL  -  LENGTH  OF  MONOPOLE  IN  WAVELENGTHS. 

C  FMC  -  FREQUENCY  IN  MEGAHERTZ. 

C  BAR  -  RADIUS  RATIO  FOR  COAXIAL  TEED  CABLE. 

C  BL  «  RADIOS  OF  CIRCULAR  DISK  IN  WAVELENGTHS  -  EPSLN/TP 

C  NSD  -  NUMBER  OF  SEGMENTS  ON  THE  DISK. 

C  NSW  -  NUMBER  OF  SEGMENTS  ON  THE  WIRE. 

C  HDL  -  HEIGHT  OF  DISK  ABOVE  THE  EARTH  IN  WAVELENGTHS. 

C  ER  -  RELATIVE  PERMITTIVITY  OT  EARTH. 

C  SIG  -  CONDUCTIVITY  OF  EARTH,  MHO/M. 

C  SET  IVCJ-1  TO  WRITE  THE  CURRENT  DISTRIBUTION  CJ(N), 

C  OR  IVCJ-0  TO  SUPPRESS  VRITE0UT  OF  CJ(N). 

C  SET  HDL  -  NEGATIVE,  FOR  MONOPOLE-DISK  IN  FREE  SPACE, 

C  OR  HDL  -  POS.  FOR  FREE  SPACE  +  FLAT  EARTH. 

C  SET  DTHD  -  NEGATIVE,  TO  SKIP  THE  GAIN  CALCULATIONS. 

C 

C  TL  -  l.D-5  FOR  EPSLN  GREATER  THAN  OR  EQUAL  0.25, 

C  -  AL.D-4  FOR  EPSLN  LESS  THAN  0.25. 

AL-.003 

BAR-3. 

BL— .  1 
ER— 4  . 

FMC— 300 . 

HDL— . 0 
HL— . 229 
SIG— . 001 
TL— AL/100 . 

IWCJ-1 

IWZ-0 
NPH-6 
NSD— 10 
NSW— 4 

NEQ-NSD+NSW-1 

AK— TP*AL 

BK— TP*BL 

HK-TP*HL 

HDK— TP*HDL 

TK— TP*TL 

OMBO— TP*FMC*1 .16 

EC— DCMPLX  (ER,  -SIG/  (OMEG*E0 ) ) 

DKD— (BK-AK) /WSD 

DKW-HX/NSW 

RH2-AK+DKD 

IT  (RH2.LT.  BAR*  AK)  OO  TO  400 
TDKD-2  .  *DKD 
CDKD-DCOS  (DSD) 

SDKD-DSIN (DSD) 

CDK-DCOS  (DRW) 

SDK— DSIW  (DKW) 

MAX— NSW- 1 
HA— WSD+1 

CALL  GFW4  (AK,  DKD,  DKW,  CDKD,  SDSD,  CDK,  SDK, TK,  IWI, NPH,  Ill) 

SIJ(1,1)-Z11 
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RIOMOND3.2 


IF (MSD.LE.l)GO  TO  100 

51- AK 

DO  60  J-2,MSD 

52- S1+DKD 

53- S1+TDKD 
Tl-AK 

DO  50  1-2, J 
T2— Tl+DKD 
T3— Tl+TDKD 

CALL  QDD (CDKD, SDKD, 81, S3,T1, T3, TK, XWZ,NPH, £22) 

SIJ (I, J)-Z22 
50  Tl— Tl+DKD 

CALL  QDM (AK,DKD, DKW,  CDKD,  SDKD,  SDK,  81,  S3, TK,  XWZ,MPH,  Z12) 

Z1 J (1, J) -Z12 
60  SI— Sl+DXD 
100  IF (MSN . LK . 1 ) GO  TO  200 

CALL  SPART  (AK,DKD,DKW,MAX,  IWZ,  ZJ,  CJ) 

L-0 

DO  160  I-NA,NEQ 
DO  150  J— I , MEQ 
K-J-X+l 

150  ZIJ(I,J)-ZJ(K) 

L— L+l 

ZZJ (1, I)-CJ(L) 

160  CONTINUE 

178  IF (NSD . LE . 1 ) GO  TO  200 
Z2-.0 

DO  190  J-NA,NEQ 
Z2-Z2+DKW 
31-22-DKW 
S3— Z2+DKW 
RH2-AK 

DO  180  I— 2 , USD 
RH2-RH2+DKD 
T1-RH2-DKD 
T3— RH2+DKD 

CALL  SXENT (AK, SI , S3, Tl, T3, CDK, SDK, CDKD, SDKD, ZKZ, Z12) 

180  ZIJ (I, J)— Z12 
190  CONTINUE 

200  CALL  GRILL  (AK,  BAR,DKD,DKW,XEQ,N3D,MSN,  TK,VJ) 

DO  210  I— 1 , NEQ 
DO  210  J-I,NEQ 
WRITE  (17,1)  I,  J,  ZIJ  (I,  J) 

ZIJ  (J,  I) -ZIJ  (I,  J) 

210  VIJ(I,J)-ZIJ(I,J) 

WRITE<17,5) 

CALL  CMINV  (CJ,VJ,  ZIJ,  ICC,  ZWCJ,  1,LLL,MM,WSQ,DET) 

Til— CJ  (1) 

Zll-l./Tll 

WRITE(6,2)AL,BL,BL,Z11 
WRITE (17 , 2 ) AL, BL, HL, Zll 

WRITE (17, 5)  _  _  . 

C  FOR  MOHOPOLE  ON  DISK  IN  FREE  SPACE,  SKIP  TO  STATEMENT  320. 
C 

r»T.T.  DZ11(AK,BAR,DKD,DKW,EC,FB,EDK,TK,IFB,D11,DV1) 

SIJ (1,1) -Dll 
IF (W8D.LE. 1) 60  TO  265 
82— AK+DKD 
DO  260  J— 2,NSD 
T2-AK+DKD 
112-1 

DO  250  X-2, J 

CALL  DSDD(AK,DBET,DKD,DKW,EC,FB,HDK, 

2, ZFB, Z12,XMX,D12,D22) 

IF (I .EQ.2) P12-D12 
ZIJ (I, J)-D22 
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112-2 

250  T2-T2+DKD 

ZZJ(1, J)-F12 
260  82-S2+DKD 
265  ZF(NSN.LE.l)GO  TO  278 
DO  276  K-1,KAX 

CALL  DZNN (AX,  DKD,  DKN,  EC,  HDK,  K,  TK, ZJ,DZ1J) 

J-NA+K-1 

ZZJ  <1,  J)  — DZ1J 

L— 1 

DO  270  I -HA, J 
ZZJ  (Z,  J)  — ZJ  (L) 

270  L-L+l 
276  CONTINUE 

278  ZF(NSD.LE.l)GO  TO  300 
Z2-.0 

DO  290  J-NA, NEQ 
Z2-Z2+DKN 

CALL  DZND  (AT, DKD, DKN, SC, HDK, NSD, TK,  Z2 ,  Z J) 

DO  280  Z-2,NSD 
280  ZZJ(Z,J)-ZJ(Z) 

290  CONTI NOT 
300  DO  310  Z— 1 , NSQ 
DO  308  J— Z , NSQ 
Z12— VXJ(Z,  J) 

D12-ZIJ(I,  J) 

MUTE  (17,1)  Z,  J,  Z12, D12 
ZZJ  (I,  J)-Z12+D12 
308  CONTINOT 
310  CONTINUE 

NRZTE (17,5) 

VJ  (1)  — VJ  (1)  +DV1 
DO  315  Z-1,NEQ 
DO  312  J-Z.NEQ 
312  ZZJ ( J,  I)  — ZZJ(Z,  J) 

315  CONTINUE 

CALL  CMZNV  (CJ, VJ,  ZZJ,  ZCC,  ZNCJ,  1,  LLL,MMM,  NEQ,  DET) 
Til— CJ (1) 

Nil— 1 . /Til 

NRZTE (6, 2) Zll , Nil 

NRZTE (17, 1) NSD, NSN,AL,BL,HDL 

NRZTE  (17 , 5) 

NRZTE (17,2) Zll, Nil 
320  CONTINUE 
400  CONTINUE 
500  CALL  BXZT 
END 
C 
C 
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TABLE  I 

DOUBLE  PRECISION 


NSD 

MSN 

AL 

BAR 

SR 

me 

SDL 

BL 

-810 

10 

4 

0.003 

3. 

4. 

300 

.0 

0.229 

.001 

BL 

ANTENNA 

IMPEDANCE  Zll 

antenna  IMPEDANCE  Ell 

(in  frmm  mpm.c%) 

(on  flat 

•axth) 

Rll 

XI 1 

Rll 

XI 1 

0.1 

14.0427 

-53.8714 

29.5665 

-32.S406 

0.2 

17.4560 

-21.1708 

25.4167 

-19.7828 

0.3 

20.3451 

-9.0883 

23.5883 

-9.5982 

0.4 

26.0578 

-1.0801 

29.5489 

-4.8751 
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TABLE  II 

DOUBLE  PRECISION 


5 

Currant  Distribution  on  aonopole  on 
Circular  Disk  in  Free  Space. 


X 

CJ(I) 

CJ(I) 

CJ(X) 

(non) 

(nag. ) 

(phase) 

1 

1.000 

0.0178959 

74.6 

2 

0.983 

0.0175990 

-105.5 

3 

0.946 

0.0169306 

-105.9 

4 

0.908 

0.0162422 

-106.1 

5 

0.857 

0.0153451 

-106.3 

6 

0.800 

0.0143135 

-106.4 

7 

0.729 

0.0130419 

-106.6 

8 

0.642 

0.0114806 

-106.7 

9 

0.528 

0.0094462 

-106.8 

10 

0.395 

0.0070747 

-106.9 

11 

0.843 

0.0150935 

72.6 

12 

0.648 

0.0116011 

71.4 

13 

0.386 

0.0069004 

70.4 

Currant  Distribution  on  Monopola  on 
Circular  Disk  on  Flat  Earth. 


X 

CJ(X) 

CJ(X) 

CJ(I) 

(non) 

(■■g.) 

(phase) 

1 

1.000 

0.0227445 

47.7 

2 

0.987 

0.0224556 

-132.2 

3 

0.957 

0.0217627 

-132.5 

4 

0.909 

0.0206744 

-132.6 

5 

0.845 

0.0192282 

-132.7 

6 

0.782 

0.0177918 

-132.6 

7 

0.713 

0.0162147 

-132.1 

8 

0.616 

0.0140152 

-131.3 

9 

0.470 

0.0107005 

-130.3 

10 

0.293 

0.0066696 

-129.4 

11 

0.880 

0.0200169 

43.9 

12 

0.686 

0.0155928 

42.0 

13 

0.410 

0.0093323 

40.4 

NSD  MSN  AL  B L  HDL 

10  4  0.3000B-02  0.1000E+00  O.OOOOE+OO 


Antenna  Xapedance  111 
(in  fraa  spaca) 

A  Z 

14.84274593  -53.87146810 


Antanna  l^adanca  Ell 
(on  Flat  earth) 

A  X 

29.58653347  -32.54060135 
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TABLE  III 
SINGLE  PRECISION 


NSD 

NS* 

XL 

BAR  ER 

PMC 

BDL 

HL 

8IG 

10 

4 

.003 

3.  4. 

300. 

.0 

0.229 

.001 

BL 

ANTENNA  IMPEDANCE  Ell 

ANTENNA  IMPEDANCE  Ell 

(in  tr¬ 

apaca) 

(on  flat  aartb) 

ull 

*11 

Rll 

*11 

0.1 

14.0184 

>53.8642 

29.5429 

-32.5342 

0.2 

17.4301 

-21.1669 

25.3907 

-19.7795 

0.3 

20.5567 

-9.0920 

23.5997 

-9.6021 

0.4 

26.8575 

-1.0805 

29.5487 

-4.6758 
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SUBROUTINE  BES10  (XX,B,B1,  ID)  BES10 

B  -  BESSEL  FUNCTION  J  sub  0  with  rwal  trguatflt  XX. 

B1  -  BESSEL  FUNCTION  J  rob  1  with  rwal  arguaant  XX. 

SET  ID  -  (0,  1,  2)  TO  CALCULATE  (J  sub  o,  J  tub  1,  or  both)  . 
IMPLICIT  REAL *8  (A-H)  ,  (P-Z) 

B-l . 

B1-.0 

IF (XX. EQ. .0) RETURN 
X-DABS (XX) 

IF (X . ST . .01) GO  TO  10 

X2-X*X 

X4-X2*X2 

B  -1.-X2/4. 4X4/64. 

Bl-X*(l.-X2/8.)/2. 

RETURN 
10  DX-X 

IF (X.GE. 3. ) GO  TO  100 
C-DX«DX/9. 

IF (ID.EQ. l)GO  TO  20 

B  -( ( ( ( (.21D-3*C- .394440-2) *C4. 444479D-1) *C-. 3163866) *C4.1265 
1 620804-1 )  *C-2 .2499997)  *C41 . 

IF (ID . EQ . 0 ) RETURN 

20  B1  -(<((((•  H09D-4*C- . 31761D-3) *C4 . 4433190-2) *C- .39542890-1) *C4 
1.21093573) »C- .56249985) *C4.5) *DX 
RETURN 
100  D»3./DX 

C-l .  /DSQRT  (DX) 

IF (ID.EQ.l)GO  TO  120 

EA-C*  ((((((- 14476D-3*D- .728050-3) *04 . 137237D-2) *D- . 9512D-4) *D- 
4 . 5527400-2) *D- .77D-6) *04.797884560803) 

FA-( ( ( ( (.135580-3*0-. 29333D-3) *0- .54125D-3) *04. 2625730-2) *D- 
5 . 3954D-4) *D- . 41663970-1) *D- .78539B1633974DX 
B  -EA*DCOS  (FA) 

IF ( ID. EQ.O) RETURN 

120  EB-C* ( ( ( ( ( (-.200330-3*04. 1136530-2) *0-. 249511D-2) *04. 17105D-3)*D4 
6.1659667D-1)*D4. 156D-5) *04.797884560803) 

FB-  ( ( ( ( (- . 29166D-3*D4 . 798240-3) *D4 . 74348D-3) *D- . 637879D-2) *04 
7 . 5650D-4) *D4 . 12499612) *0-2 . 3561944901924DX 
B1  -EB*DC0S(FB) 

RETURN 

END 
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SUBROUTINE  CISI <CI,CIN, SI,X)  CISI 

CALCULATES  CZ  -  COSINE  INTEGRAL,  AND 
SI  -  SINE  INTEGRAL  WITH  ARGUMENT  X. 

IMPLICIT  REAL'S  (A-B) ,  (P-Z) 

DATA  GAM.P2/. 57721566, 1.57079632/ 

A-DABS (X) 

IF (A.GT . 4 . ) GO  TO  10 

IF (A.GT. . 1) GO  TO  3 

IF(A.GT.O.)GO  TO  2 

CI-.O 

CIN- . 0 

SI-.O 

RETURN 

2  X2-A*A 

SI-X* ( ( . 03*X2-1 . ) *X2/18 .+1 . ) 

CIN- . 25*X2* ( (X2/45 . -1 . ) »X2/24 . +1 . ) 

GO  TO  B 

3  Y-(4.-A)*(4.+A) 

SI-X*  < ( ( ( (1 . 753141D-9*Y+1 . 568988D-7) *1+1 . 374168D-5) *Y+6 . 939889D-4) 
C*  Y+l . 964882D-2) *Y+4 . 395509D-1) 

CIN-  A*A* ( < ( (  (1 . 386985D-10*Y+1 .584996D-8) *T 

C+l -725752D-6) *Y+1 . 185999D-4) *Y+4 . 990920D-3) *Y+1 . 315308D-1) 

8  Cl— GAM+DLOG (A) -CIN 

RETURN 

10  SI-DSIN(A) 

Y-DCOS (A) 

Z— 4 . /A 

O-  ((((((( <4 . 048069D-3**-2 -279143D-2) *Z+S .515070D-2) *Z-7 . 261642D-2) 
C*Z+4 . 987716D-2) *Z-3 . 332519D-3) *Z-2 . 31461 7D-2) *Z-1 . 13495BD-5) *Z 
C+6.250011D-2) 'Z+2.583989D-10 
V-( (((((( ( (-5.108699D-3*Z+2.819179D-2)  *Z-6.537283D-2) *Z 
C+7 . 902034D-2) *Z-4 . 400416D-2) *Z-7 . 94SS56D-3) *Z+2 . 601293D-2) *Z 
C-3.764000D-4) *Z-3 . 12241BD-2) *Z-6.646441D-7) *Z+2.5D-1 
CI-Z* (SI*V-Y*U) 

SI— Z*  (SI*U+Y*V)+P2 
IF  (X.LT.  .0)SI— SI 
CIN-GAM+DLOG (A) -Cl 
RETURN 
END 
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CMINV.  1 


e 

SUBROUTINE  CMINV  (C,V,  Z,  IDM,  INR,  112 ,  L, M, NEQ, DET) 

C  CMINV2  INVERTS  THE  MATRIX  Z  (I,  J)  AND  SOLVES  THE 

C  SIMULTANEOUS  LINEAR  EQUATIONS  TO  DETERMINE  C(I) . 

C  V(I)  -  EXCITATION  COLUMN. 

C  Z(I,J)  -  IMPEDANCE  MATRIX. 

C  IDM  -  DIMENSION  07  Z  (IDM,  IDM)  IN  CALLING  PROGRAM. 

C  INR  >  1  IF  SOLUTION  IS  TO  BE  PRINTED. 

C  INR  >  0  H  PRINTOUT  IS  TO  BE  SUPPRESSED. 

C  112  -  1  ON  FIRST  CALL,  NKERE  CMINV  MUST  INVERT  Z. 

C  112  ■  2  ON  LATER  CALLS,  IF  Z(I,J)  BAS  ALREADY  BEEN  INVERTED. 

C  L(I) ,  M(I)  -  NORX  ARRAYS. 

C  NEQ  -  NUMBER  OF  SIMULTANEOUS  LINEAR  EQUATIONS . 

C  DET  -  DETERMINANT  OF  TBE  SQUARE  MATRIX. 

IMPLICIT  REAL* 8  <A-H) ,  (P-Z) 

COMPLEX* 16  C(1),V(1),S 

COMPLEX* 16  Z (IDM, IDM), BIGZ, BOLD, DET 

DIMENSION  L(l) , M fl) 

2  rORMAT (IX, I5,F10 .3,715.7,710.1) 

5  FORMAT (1B0) 

N-NEQ 

IF (112 .NE . 1) GO  TO  150 
C  DET-DCMPLX (1.00,0. DO) 

DO  80  K**1,N 
L(K)-K 
M(K)-K 
BIGZ-Z (K,K) 

DO  20  J-K,N 
DO  20  I-K,N 

10  IF (CDABS (BIGZ) -CDABS  (Z (I, J) ) ) 15,19,19 
15  BIGZ-Z  (I,  J) 

L  (K)  —I 
M(K)-J 

19  CONTINUE 

20  CONTINUE 
J-L(K) 

IF (J-K) 35,35,25 
25  CONTINUE 

DO  30  I— 1,N 
BOLD—  Z(IC,I) 

Z  (K,  I)  — Z  ( J,  I) 

30  Z ( J, I) —BOLD 
35  I-M  (K) 

IF(I-K)  45,45,38 
38  CONTINUE 

DO  40  J-1,N 
BOLD— Z(J,K) 

Z(J,K)  -Z(J,I) 

40  Z (J, I) —BOLD 
45  CONTINUE 

DO  55  I— 1,N 
IF (I-K) 50,55,50 
50  Z  (I,K)— Z  (I,K)  /  (-BIGZ) 

55  CONTINUE 

DO  65  I— 1,N 
DO  85  J-1,N 
IF  (I-K)  60,  64,  60 
60  IF  (J-K)  62,  64,  62 
62  Z  (I,  J)  — Z  (X,K)  *Z  (K,  J)  +Z  (I,  J) 

64  CONTINUE 

65  CONTINUE 

DO  75  J— 1,N 
IF (J-K) 70,75,70 
70  Z  (K,  J)  — Z  (K,  J)  /BIGZ 
75  CONTINUE 
C  DET— DET*BIGZ 
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Z(K,K)-1./BIGZ 
80  CONTINUE 
K-N 

100  K-K-l 

ir (K) 150, 150, 105 
105  I-L(K) 

IF(I-K) 120,120,108 
108  CONTINUE 

DO  110  J-1,N 
HOLD-Z(J,K) 

Z(J,K)—  Z(J,I) 

110  Z  ( J,  I) “HOLD 
120  J-M(K) 

IF (J-K) 100, 100, 125 
125  CONTINUE 

DO  130  I-1,N 
HOLD-Z (K, I) 

Z(K,D— Z(J,I) 

130  Z(J,I)-HOLD 
GO  TO  100 
150  CMX-.O 

DO  220  I-1,NEQ 
S-DCMPLX ( . 0D0 , .ODO) 

DO  210  J-1,NEQ 
210  S-S+Z(I, 

SA-CDABS (S) 

IF (SA . GT . CKX) CHX-SA 
220  C(I)-S 

IF ( IWR . LE . 0 ) GO  TO  250 
WRITE  (17 , 5) 

DO  240  I-1,NEQ 
S-C(I) 

SA-CDABS (S) 

SN-SA/CMX 

PH-.O 

IF(SA.LE. .0)GO  TO  240 
PH-57. 29578 *DATAN2  (DIMAG (S)  , DREAL (S)  ) 
240  WRITE(17,2)I,SN,SA,PH 
WRITE (17,5) 

250  RETURN 
END 
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SUBRODTIME  DZ11  (AK, BAR, DKD, DKN, EC, FB, HDK,  TK,  ZFB,D11,DVX) 

C  DZ11  CALCULATES  Dll  -  CHANGE  IN  SELF- IMPEDANCE  OF  MODE  1 
C  DUE  TO  REFLECTION  FROM  FLAT  EARTH. 

C  ALSO  DV1  -  ONE  TERM  IN  VOLTAGE  FOR  MODE  1. 

IMPLICIT  REAL*8  (A-H) ,  (P-Z) 

DIMENSION  FB (1) 

C0MPLEX*16  FST,G,GAM,RC,EC,ZAA,ZHH,PC,EG1,EG2 
COMPLEX* 16  ZDD, ZMM, EST, ERD, ERM, EGZ, CB, Dll , ZDM, ZMD 
COMPLEX* 16  DV1, VDD, VDM,VHH, VAA 

DATA  ETA, PI, TP/376. 730366239, 3. 14159265359, 6. 28318530718/ 
BAL-DLOG (BAR) 

BK-BAR*AK 
RH2-AK+DKD 
CDK-DCOS (DKN) 

SDK-DSIN (DKW) 

SDKD-DSIN (DKD) 

DK1-DKW+TK 

SDKl-DSIN(DKl) 

CDK1-DCOS (DK1) 

EST-DCMPLX ( . 0D0, -ETA/ (4 . *PI*SDK1) ) 

FST-DCMPLX ( . 0D0, -ETA/  (4 . *PI*SDKD) ) 

DBET-.l 

KMX-200 

IF (KMX.GT.IFB)KMX-IFB 

C  NEXT  CALCULATE  F (BETA)  BY  INTEGRATING  ACROSS  THE  DISK. 

DO  60  K— 1 , KMX 
DRK-PI/10. 

BET-DBET* (K-l) 

IF (BET . GT . 1 . ) DRK-DRK/BET 
I NT- (RH2-AK) /DRK 
IF (XNT .LT . 10) INT— 10 
DRK- (RH2-AK) /I NT 
F-.O 

RK— AK+DRK/ 2 . 

C  NEXT  INTEGRATE  ACROSS  THE  DISK. 

DO  50  1-1, INT 
BRhBET * RK 

CALL  BES10  (BR, BJ0 ,BJ1, 0) 

F-F+BJ0*DCOS (RH2-RK) 

50  RK-RK+DRK 
FB (K) — DRK*F 
60  CONTINUE 
C  NEXT  CALCULATE  ZDD. 

ZK-HDK+TK 
Zl— HDK 

Z2-HDK+TK+DKW 
ERD-DCMPLX ( . 0D0, . ODO) 

ERM-DCMPLX ( . ODO , . ODO ) 

VDD-DCMPLX ( . ODO , . ODO) 

C  NEXT  INTEGRATE  ON  BETA. 

DO  80  K-l, KMX 
BET-DBET* (K-l) 

F-FB (K) 

CALL  BES10  (BET*AK,BA0,BA1, 0) 

BETS— BET  *BET 

IF (BET . GT . 1 . ) GO  TO  62 

BR-DSQRT (1 . -BETS) 

ARG— HR*ZK 

EGZ— D  CMP  LX  (DCOS  (ARG) ,  -DSIN  (ARG) ) 

ARG— HR*Z1 

EG1— DCMPLX  (DCOS  (ARG)  ,  -DSIN  (ARG) ) 

ARG— HR* Z 2 

EG2— DCMPLX  (DCOS  (ARG) ,  -DSIN  (ARG)  ) 

GAM— DCMPLX ( . ODO , HR) 
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GO  TO  64 

62  ALP-DSQRT (BBTS-1 . ) 

EG2-DCMPLX  ( . ODO, . 0D0) 

ARG-ALP*Z2 

IF (ARG .LT. 80 . ) EG2 -DCMPLX (DEXP (-ARG) , .0D0) 
ARG-ALP*ZK 

IF (ARG.GT.BO.)GO  TO  81 
EGZ-DCMPLX (DEXP (-ARC) , . 0D0) 

EG1-DCMPLX (DEXP (-ALP*Z1) , . 0D0) 

GAM -DCMPLX  (ALP,  .  0D0) 

64  G— CDSQRT  (BETS -EC) 

RC- (GAM*EC-G) / (GAM*EC+G) 

CB— EG2 -CDK1 *EG1 
C  NEXT  INTEGRATE  ON  RRO. 

RSJ- . 0 
DRX-PI/10. 

IF (BET . GT . 1 . ) DRK-DRK/BET 

INT-DKD/DRK 

IF (INT . LT . 10) INT-10 

DRK-DKD/INT 

RK-AK+DRK/2 . 

DO  70  I-l, INT 

CALL  BES10  (BET*RK, BJO , BJ1 , 1) 

RSJ-RSJ+BJ1 *DS  IN (RH2-RX) 

70  RK-RK+DRK 
RSJ-DRK*RSJ 

ERD-ERD+RS J*GAM*RC*F»EG1 *EGZ 
ERM-ERM+RS  J*RC*BA0  *CB*KGZ 
CALL  BES10  (BET*BK,BB0,BB1, 0) 

VDD-VDD+RSJ»RC* (BAO-BBO) *EGZ*EGZ 

80  CONTINUE 

81  ERD-DBET*ERD*DCMPLX  ( .  ODO ,  -ETA/  (4  .  *PI*SDKD*SDKD)  ) 
ERM-DBET*ERM*DCMPLX  ( .  ODO,  -ETA/  (4 .  «PI*SDKD*SDK)  ) 
ZDD-ERD+ERM 

VDD— DBET*VDD/ (2 . *BAL*SDKD) 

C  NEXT  CALCULATE  ZDM  BY  INTEGRATING  ON  BETA. 
Zl-HDK+TK 
Z2-Z1+DKN 

ZDM-DCMPLX ( . ODO , . ODO) 

DO  90  K-l.KMX 
BET-DBET* (K-l) 

F-FB (K) 

BA-BET* AT 

CALL  BES10  (BA, BJO, BJ1 , 0) 

BETS-BET*BET 

IF (BET . GT . 1 . ) GO  TO  82 

HR-D30RT (1 . -BETS) 

ARG^HR^HDK 

EGZ-DCMPLX (DCOS (ARC) , -DSIN (ARC) ) 

ARC— BR*Z1 

EG1— DCMPLX (DCOS (ARG) , -DSIN (ARG) ) 

ARG— HR* Z 2 

BG2— DCMPLX (DCOS (ARG) , -DSIN (ARG) ) 

GAM— DCMPLX ( . ODO , HR) 

•  GO  TO  84 

82  ALP-DSQRT (BETS -1.) 

EG2— DCMPLX ( . ODO , . ODO) 

ARG— ALP* Z2 

IF  (ARC .  LT .  80 . )  EG2-DCMPLX  (DEXP  (-ARG) ,  .ODO) 

ARG— ALP*Z1 

IF (ARG. GT. 80.) GO  TO  92 
EG1— DCMPLX (DEXP (-ARC) , .ODO) 

EGZ-DCMPLX (DEXP (-ALP*BDK) , .ODO) 

GAM— DCMPLX (ALP , .ODO) 

84  G— CD3QRT (BETS-EC) 

RC- (GAM*EC-G) / (GAM*EC+G) 
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ZF (K.GT.l)GO  TO  86 

PC-  (DCMPLX  (1  .D0,2  .DO*DKW)  *EG2-DCMPLX  (CDK,  SDK)  *EG1> /A  . 
GO  TO  90 

86  PC«((GAM*SDK-CDK)*EG1+EG2)/BETS 
90  ZDM«ZDM+BET*RC*F*BJO*EGZ*PC 
92  ZDM— DBET*FST*ZDM/SDK 
C  NEXT  CALCULATE  ZHH  BY  INTEGRATING  ON  H. 

VHH— DCMPLX ( . 0D0 , . 0D0 ) 

ZHH— DCMPLX ( . 0D0 , . ODO ) 

DH-DBET 
NS— 1 . /OH 
DH-1 . /NS 
HR-DH/2 . 

DO  100  I— 1,NS 
BETS— 1 . -HR*HR 
BET-DSQRT (BETS) 

CALL  BES10  (BET*AK,BA0,BA1,0) 

CALL  BES10  (BET*BK,BBO,BB1,0) 

GAM— DCMPLX  (  .  ODO  ,  HR) 

G— CDSQRT (BETS-EC) 

RC- (GAM*EC-G) / (GAM*EC+G) 

EGZ— DCMPLX  (DCOS  (HR*HDK)  ,  -DSIN  (HR*HDK)  ) 

EG1— DCMPLX (DCOS  (HR*Z1) , -DSIN(HR*Z1) ) 

EG2— DCMPLX (DCOS  (HR*Z2) , -DSIN (HR»Z2)  ) 

PC- ( (GAM* SDK- CDK) *EG1+EG2) /BETS 
HR— HR+DH 

VHH— VHH+RC*  (BAO-BBO)  *BA0*EG1*PC 
100  ZHH-ZHH+RC*BAO*BAO*  (EG2-CDK1*EGZ)  *PC 

VHH— DH*VHH*DCMPLX  ( .  ODO ,  -1  .DO/  (2  .D0*BAL*SDK)  ) 

ZHH— DH*ZHH*ETA/  <4  .  »PI»SDK»SDK1) 

C  NEXT  CALCULATE  ZAA  BY  INTEGRATING  ON  ALP. 

DA-DBET 
ALP-DA/2 . 

VAA— DCMPLX ( . ODO , . ODO ) 

ZAA— DCMPLX ( . ODO, . ODO) 

DO  110  I— 1, KMX 
BETS— ALP* ALP +1 . 

BET-DSQRT (BETS) 

CALL  BES10  (BET*AK,BA0,BA1,0) 

CALL  BES10  (BET*BK, BBO , BB1 , 0) 

G— CDSQRT (BETS-EC) 

RC- (ALP*EC-G) / (ALP*EC+G) 

KA2-.0 
ARG— ALP*Z2 

IF  (ARG . LT .  80 . )  EA2-DEXP  (-ARG) 

ARG— ALP*Z1 

IF (ARG . GT . 80 . ) GO  TO  112 
EA1-DEXP (-ARC) 

EAZ-DEXP (-ALP*HDK) 

P-  ( (ALP  *  SDK-  CDK )  *EA1+EA2)  /BETS 
VAA— VAA+RC*  (BAO-BBO)  *BA0*EA1*P 
ZAA— ZAA+RC*BA0  *BA0  *  (EA2-CDK1*EAZ)  *P 
110  ALP-ALP +DA 
112  ZAA— DA*EST  *  ZAA/ SDK 

VAA— DA*VAA/ (2 . *BAL*SDK) 

ZMD— ZAA-f  ZHH 
ZMM— ZDM-fZMD 
DV1-VAA+VDD+VHH 
Dll-ZDD+ZMM 
RETURN 
END 
C 
C 
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c 


c 


c 

c 


SUBROUTINE  DZDD (AK, DBET, DKD,DKW,EC,FB, HDK, S2, T2,TK  DZDD.l 

2, IFB, 112, KMX, D12, D22) 

DZDD  CALCULATES  D12  -  CHANGE  IN  MUTUAL  IMPEDANCE  BETWEEN  MODE  1 
AND  DISK  DIPOLE  MODE. 

ALSO  D22  -  CHANGE  IN  MUTUAL  IMPEDANCE  BETWEEN  TWO  DISK  DIPOLE  MODES. 
IMPLICIT  REAL'S  (A-H) ,  (P-Z) 

DIMENSION  FB (1) 

COMPLEX* 16  D12,D22,EC,QST,GAM,G,RC,EGZ,ZDD,ZDW,EZD,EZ1,EZ2,PC,QDW 

DATA  £~A,PI, TP/376. 730366239, 3. 14159265359, 6. 28318530718/ 

SDKD-DSIN (DKD) 

CDKD-DCOS (DKD) 

S1-S2-DKD 
S3— S2+DKD 
T1-T2-DKD 
T3-T2+DKD 

QST-DCMPLX  ( .  0D0 ,  -ETA /  (4 . *PI*SDKD*SDKD) ) 

IF(I12.GT.l)GO  TO  62 

DBET-.l 

KMX-200 

IF (KMX . GT . IFB) KHX-IFB 

NEXT  CALCULATE  F  (BETA)  BY  INTEGRATING  ACROSS  THE  ANNULAR  DISK. 

DO  60  K— 1, KMX 
DRK— PI/10 . 

BET-DBET* (K-l) 

IF (BET . GT . 1 . ) DRK— DRK/BET 
INT— DKD/DRK 
IF (INT .LT . 10) INT— 10 
DRK— DKD  /  INT 


F-.O 

RK— Sl+DRK/2 . 

NEXT  INTEGRATE  ACROSS  THE  DISK. 

DO  50  L-1,2 
DO  40  I— If INT 

CALL  BES10  (BET'RX, B JO , B J1 , 0) 

IF (L . EQ . 1 ) F-F-B JO  *DCOS (RK-S1 ) 

IF (L . EQ . 2 ) F-F+BJO*DCOS ( S3-RK) 

40  RK-RK+DRK 
50  RK— S2+DRK/2 . 

FB (K) — DRK*F 
60  CONTINUE 
62  CONTINUE 
NEXT  CALCULATE  D22 . 

ZK— 2 . *HDK+TK 
NEXT  INTEGRATE  ON  BETA. 

D22— DCMPLX(.0D0, .000) 

DO  100  K-l, KMX 
DRK— PI/10 . 

BET-DBET* (K-l) 

IF (BET . GT . 1 . ) DRK— DRK/BET 
INT— DKD /DRK 
IF  ( INT .  LT .  10 )  INT-1 0 
DRK— DKD/  INT 
F— FB (K) 

BETS— BET'BET 

IF (BET . GT . 1 . ) GO  TO  72 

HR-DSQRT (1 . -BETS) 

ARC— HR'ZK 

EGZ— DCMPLX  (DCOS  (ARC)  ,  -OSIN  (ARC)  ) 
GAM— DCMPLX ( . 0D0 , HR) 

GO  TO  74 

72  ALP— DSQRT (BKTS-1 . ) 
arg-alp*zk 

IF (ARG . GT . 80 . ) GO  TO  102 
EGZ-DCMPLX (DEXP (-ARG) , .0D0) 
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GAM-DCMPLX  (ALP ,  .  ODO) 

74  G-CDSQRT (BETS-EC) 

RC- (GAM*EC-G) / (GAM*EC+G) 

C  NEXT  INTEGRATE  ACROSS  THE  ANNULAR  DISK. 
R22-.0 
RK-T1+DRK/2. 

DO  90  L-1,2 
DO  80  I— 1 , INT 

CALL  BES10  (BET*RX,BJ0,BJ1,1) 
ir (L .EQ . 1 >  R22-R22+BJ1*DSIN  (RK-T1) 
ir (L. EQ . 2) R22-R22+BJ1 *DSIN (T3-RK) 

80  RK-RK+DRK 
90  RK-T2+DRK/2 . 

100  D22-D22+DRK*GAM*RC*F*EGZ*R22 
102  D22-DBET*QST*D22 

IF (112 .NE . 1) RETURN 
c  NEXT  CALCULATE  D12. 

ZDD-DCMPLX ( . ODO, . ODO) 

ZDN-DCMPLX ( . ODO , . ODO) 

R2-AK+DKD 

C  NEXT  INTEGRATE  ON  BETA. 

Zl-HDK+TK 
Z2-Z1+DKN 
SDK-DSIN (DKH) 

CDK-DCOS  (D1CN) 

QDH-DCMPLX  ( .  ODO,  -ETA/  (4 .  *PI*SDX*SDKD)  ) 

DO  160  K-1,KMX 

DRX-PI/10. 

BET-DBET* (K-l) 

IT  (BET .  GT .  1 .  )  DRK-DRK/BET 

INT-DKD/DRK 

IF  (INT  .LT .  10)  INT— 10 

DRK-DKD/INT 


F-EB  (K) 

BETS«BET*BET 

IF (BET . GT . 1 . ) GO  TO  112 

HR-DSQRT (1 . -BETS) 

ARG-HR*ZK 

EGZ-DCHPLX  (DCOS  (ARC)  ,  -DSTN  (ARC)  ) 
ARG-HR*HDK 

EZD-DCMFLX  (DCOS  (ARG)  ,  -DSIN  (ARG) ) 

ARG-HR*Z1 

EZ1-DCMPLX  (DCOS  (ARG)  ,  -DSIN  (ARG)  ) 

ARG-HR*Z2 

EZ2-DCMPLX  (DCOS  (ARC) ,  -DSIN  (ARG)  ) 
GAM-DCMPLX ( . 0D0 , HR) 

GO  TO  114 

112  ALP-DSQRT (BETS-1 . ) 

EZ2-DCMPLX(.0D0, .0D0) 

ARG-ALP*Z2 

IF(ARG.LT.80.)EZ2-DCMPLX(DEXP(-ARG),  .0D0) 
EGZ-DCMPLX ( . ODO , . ODO) 

ARG*AXf  *ZK 

IF  {ARG.  LT.  80.)  EGZ-DCMPLX  (DEXP  < -ARG),  .0D0) 

ARG-ALP*Z1 

IF (ARG.GT.80 . ) GO  TO  162 
EZ1-DCMPLX (DEXP (-ARC) , .ODO) 

EZD-DCMFLX  (DEXP  (-A1PUDK)  ,  .0D0) 
GAM-DCMPLX  (ALP,  .0D0) 

114  G-CDSQRT (BETS -EC) 

RC- (GAM*BC-G) / (GAM*EC+G) 

C  NEXT  INTEGRATE  ACROSS  THE  DZSX. 

R22-.0 

RK-AK+DRX/2 . 

DO  140  I— 1, INT 

CALL  EES 10  '  (BET*RK,BJ0, BJ1, 1) 
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R22-R22+BJ1 *DSIN (RX-R2) 

140  RK-RK+DRK 

CALL  BBS10  (BET*AK, BJO , BJ1 , 0 ) 

IF<K.GT.l)GO  SO  150 
PC-  (DOS’ LX  (1  .DO,  2 .  *DKW)  *ZZ2-DCMPLX  (CDK,  SDK)  *KZ1)  /4  . 
GO  TO  152 

150  SC— ( (GAM*  SDK-CDK) *KZ1+EZ2 ) /BETS 
152  ZDM— ZDW+BET*RC*F*BJO*EZD*PC 
160  ZDD-ZDD+DRK*GAM*RC*F*EGZ*R22 
162  ZDD— DBET*QST*ZDD 
ZDW-DBET * QDB *  ZDM 
D12— ZDD+ZDW 
RETURN 
END 


A-20 


c 


DZWD.l 


16 

SUBROUTINE  DZWD  <AK, DKD, DKW,EC, HDK,NSD,  TK,  ZK2,  DZ12) 

C  DZWD  CALCULATES  DZ12  -  CHANGE  IN  MUTUAL  IMPEDANCE 

C  BETWEEN  A  WIRE  DIPOLE  MODE  AND  A  DISK  DIPOLE  MODE. 
IMPLICIT  REAL*6  (A-H)  ,  (P-Z) 

COMPLEX*  16  DZ12<1) ,  EC, GAM, G, RC, QC,BGZ, QST 

DATA  ETA,PI, TP/376. 730366239, 3. 14159265359, 6.28318530718/ 

IF (NSD.LE.l) RETURN 

DO  20  1-2, NSD 

20  DZ12 (I) — DCMPLX ( . 0D0, . 0D0) 

TDKD— 2 . *DKD 
CDK-DCOS (DKW) 

SDK— DS  IN  (DKW) 

SDKD-DSIN (DKD) 

Zl-HDK+TK 
ZJ2-Z1+ZK 2 
ZT-Z1+ZJ2 
DBET-.l 
KMX-200 
BET— DBET/2 . 

C  NEXT  INTEGRATE  ON  BETA. 

DO  X00  K— 1,KMX 
BETS— BET*BET 
IF (BET . GT . 1 . ) GO  TO  42 
BR-DSQRT (1 . -BETS) 

GAM— DCMPLX ( . ODO , HR) 

CGD-DCOS (HR* DKW) 

ARG— HR*ZT 

EGZ— DCMPLX  (DCOS  (ARG)  ,  -DSIN  (ARG) ) 

GO  TO  44 

42  ALP— DSQRT (BETS -I . ) 

GAM— DCMPLX (ALP, .0D0) 

EAD-DEXP (ALP* DKW) 

CGD— (EAD+1 . /EAD) /2 . 

ARG— ALP* ZT 

IF (ARG . GT . SO . ) GO  TO  102 
EAZ-DEXP  (-ARG) 

EGZ— DCMPLX (EAZ, .ODO) 

44  G— CDSQRT (BETS-EC) 

RC- (GAM*EC-G) / <GAM*EC+G) 

CALL  BES10  (BET*AK, BJ0 , BJ1 , 0) 

CC— CGD-CDK 
DRX-PI/10 . 

IF (BET . GT . 1 . ) DRK— DRX/BET 
ZNT— DKD/DRK 
IF (INT.LT. 10) XNT-10 
DRK-DKD/ ZNT 
RBI— AK 

QC- (RC*EGZ) * (DRK*BJ0*CC) 

DO  SO  1-2, NSD 
RH2-RH1+DKD 
RB3-RB1 +TDKD 

C  NEXT  INTEGRATE  ACROSS  THE  ANNULAR  DISK. 

FR-.O 

RX— RH1+DRK/2 . 

DO  70  L-1,2 
DO  <0  J-1,ZNT 

CALL  BSS10  (BET*RX,BJ0,BJ1,1) 

IF (L.EQ. 1)CI-DSIN (RX-RH1) 

ZF (L.EQ.2) Cl— DSIN (RH3-RK) 

FR-FR+8J1*CI 
60  RK-RX+DRK 
70  RK-RH2+DRK/2 . 

RH1-RH1+DKD 

B0  DZ12 (Z)— DZ12 (I)+QC*FR 
100  BET— BET+DBET 
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102  QST-DCMPLX  ( .  ODO ,  DBET*ETA/  (TF*SDKD*SDK) ) 
DO  120  X>2,NSD 
120  DZ12<I)-QST*DZ12(I) 

RETURN 

KND 

C 

C 


noon  on  oono 
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IS 

SUBROUTINE  DZWW (AX, DKD, DKW, EC, HDK, J, TK, DZI J,DZ1 J) 

DZMW  CALCULATES  DZIJ  -  CHANGE  IN  MUTUAL  IMPEDANCE  OF 
TWO  MIRE  DIPOLE  MODES. 

AND  DZIJ  -  CHANGE  IN  MUTUAL  IMPEDANCE  BETWEEN  MODE  1  AND 
A  WIRE  DIPOLE  MODE. 

IMPLICIT  REAL* 8  (A-H)  ,  (P-Z) 

COMPLEX* 16  DZIJ (1) , DZIJ 

COMPLEX*16  G,GAM,EC,RC,QC,QJ,EJH,EIH,EZ1,EZ2,PC,ZDW, EGZ 
DATA  ETA, PI, TP/376. 730366239, 3. 14159265359, 6. 28318530718/ 
QJ-DCMPLX ( . 0D0, 1 .DO) 

CDK-DCOS (DRW) 

SDK-DSIN (DXW) 

SDKD-DSIN (DKD) 

DO  20  I-1,J 

20  DZIJ ( I ) —DCMPLX ( . 0D0 , . 0D0 ) 

DZI  J-DCMPLX  ( . 0D0, .0D0) 

Zl-HDK+TK 
-2-Z1+DKW 
ZJ2-ZX+J»DKW 
NEXT  INTEGRATE  ON  H. 

DH-.25/ZJ2 
DH-.l 
KMX-1 . /DH 

IF (KMX. LT. 10) KMX-10 
DH— 1 . /KMX 
HR— OH/2 . 

DO  100  K-1,XMX 
BETS— 1 . -HR*HR 
BET— DSQRT (BETS) 

CALL  BES10  (BET*AK , BJ0 ,  BJ1 , 0 ) 

GAM— DCMPLX ( . 0D0 , HR) 

G— CDSQRT (BETS-EC) 

RC-  (GAM*EC-G)  /  (GAM*EC+C) 

CC-DCOS (HR*DKW) -CDK 
rAC-BJ0*CC 

QC-RC* (FAC*FAC*DH/BETS) 

ARC— HR*ZJ2 

EJH-DCMPLX  (DCOS  (ARG) ,  -D3IN  (ARC)  ) 

QC— QC*EJH 
ZZ2-Z1 
DO  80  I-l.J 
ZI2-ZI2+DKW 
ARG— HR* ZI 2 

EIH-DCMPLX  (DCOS  (ARG) ,  -DSIN  (ARG)  ) 

80  DZI J (I) —DZI J (I) +QC*EIH 
BJ2-BJ0 *BJ0 /BETS 
ARG— BR*Z1 

EZ 1— DCMPLX  (DCOS  (ARG) ,  -DSIN  (ARC) ) 

ARG— HR* Z 2 

EZ2— DCMPLX  (DCOS  (ARG) ,  -DSIN  (ARG)  ) 

PC- EZ2+EZ1  *  DCMPLX  ( -CDK,  HR*  SDK) 

DZI J— DZI J+RC*BJ2*CC*EJH*PC*DH 
100  KR— HR+DH 
NEXT  INTEGRATE  ON  ALP. 

AKX-6 . / (ZJ2+Z1) 

DA-.1/ZJ2 
KMX— AMX/DA 
DA-DH 
KMX-200 
ALP-DA/2. 

DO  200  X-1,KMX 
BET S -ALP* ALP +1 . 

BET— DSQRT (BETS) 

CALL  BES10  (BET*AK,BJ0,BJ1,0) 

G— CDSQRT (BETS-EC) 
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*C-  (XLP*EC-G)  /  (ALP*EC+G) 
ead-dsxf (ALP*DKW) 

CAD* (EAD+1 . /BAD) / 2 . 

CC-CAD-CDK 

QC^’qJ*FC?* (rAC*rAC*DA/BETS) 
XRG-ALP*ZJ2 

ir <xrg.gt.80.)go  so  202 

AJZ-DEXP (-ARG) 

qc-qc*ajz 

ZI2-Z1 

DO  180  I-lf J 
ZI2-ZX2+DKH 
XIZ-DZX3?  (-AU?*ZI2) 

180  OZIJ(I)“D2IJ<I)+QC*A2Z 
BJ2-BJ0*BJ0/B*TS 

AZ1-DEXP (-ALP*Z1) 

AZ2-DEXP ( -ALP *22) 

PR-AZ2+ (ALP*SDK-CDK) *AZ1 

DZI^-DZIJ+QJ'RC'B^’CC*^*®11  ^ 

200  ALP-ALP+DA 
202  RIJ-ETA/ (PI*SDK*SDK) 
DZlJ-RIJ*DZlJ/2 . 

DO  220  1*1/ J 
220  DZIJ(l)-RIJ*°*IJi£* 

c  HEXT  INTEGRATE  OH  BETA. 

KMX-200 

ZT-Z1+ZJ2 

RH2-AK+DKD 

xdh-dcmplx ( . odo , . ooo) 

DBET- .  1 
BET-DBET/2 . 

DO  300  K“  1,K>DC 

BETS-BET*BET 

IF (BET.GT .1 . )GO  TO  262 
HR-DSQRT (1 . -BETS) 

XGZ-DCKPLX  (DCOS  (ARG)  ,  -OSIH  (ARG) ) 
CAK-DCMPLX ( . ODO , HR) 

CGD-DCOS (HR*DKW) 

CO  TO  264 

262  ALP“DSQRT <BETS-1 . ) 

ARG-ALP*ZT 

ir(ARG.GT.BO.)GO  TO  302 
AGZ-DEXP <-ARG) 

*G2-DCHPLX{AGZ, .000) 

GJkM*DCMPLX  (ALP ,  .000) 

KAD*DEXP (ALP*DKH) 

CCD— (BAD+1 . /BAD) /2 . 

264  C-CDSQRT 

RC- (GAH*EC-G)  /  (GAM*EC+G) 

f^T.T.  BXS10  (BET*AX,BJ0,BJ1/ 0) 

OOCGD-CDK 

0ORC*BJ0*CC*EGZ 

iImBET^CT  !  1 . )  DRJC-DRX/MT 

IHT-DKD/DRX 

IF (IMS .LT . 10) IHT-10 

ORKOKD/IHT 

RK*AR+ORX/ 2 . 

RDH-.O 

C  NEXT  IHTEGRATE  OH  RHO. 

SiIkwo'^Sw^^bjI/I) 

RDH"ROH+BJl*OSm  (RH2-RK) 

280  RK-RK+ORK 
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ZDH«ZDK+DRX*QC*RDH  DZW.3 

300  BET-BET+DBET 

302  ZDN*>DBET*ZDM*DCKPIJC  ( .  0D0,  -ETA/  <TP*8BKD*SDK)  ) 

DZ1 J«DZ1 J+ZDW 

RETURN 

END 
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SUBROUTINE  EXP J (VI, V2,W12) 

C  EXP  J  CALCULATES  ¥12  -  EXPONENTIAL  INTEGRAL  WITH 

C  LOWER  LIMIT  VI  AND  UPPER  LIMIT  V2. 

IMPLICIT  REAL* 8  (A-H) ,  (P-Z) 

COMPLEX* 16  EC, E15, S, T,UC,VC, V1,V2,W12,  Z 
DIMENSION  V(21),W(21),D(16),E(16) 

DATA  V/  0 . 22284667D  00, 

20 . 11889321D  01,0.299273630  01,0.577514360  01,0.983746740  01, 

20 . 15982874D  02, 0 . 933078120-01, 0.49269174D  00,0.121559540  01, 
20.22699495D  01, 0 . 36676227D  01,0.542533660  01, 0.75659162D  01, 
20.10120228D  02,0. 13130282D  02, 0 . 16654408D  02, 0 .20776479D  02, 
20.25623894D  02, 0 . 31407519D  02, 0.38530683D  02, 0 . 48026086D  02/ 
DATA  W /  0 . 45896460D  00, 

20. 41700083D  00, 0 . 11337338D  00, 0 .10399197D-D1, 0 .26101720D-03, 

20. 8 98 5479 ID-06, 0.21823487D  00, 0 .34221017D  00, 0 .26302758D  00, 
20.12642582D  00, 0 .40206865D-01, 0 -85638778D-02, 0 . 12124361D-02, 

20. 11167440D-03, 0 . 645992670-05,0. 222631 6 9D-06, 0 . 42274304D-08, 

20 .39218973D-10, 0 .14565152D-12, 0 .14830270D-15,  0 .160059490-19/ 
DATA  D/  0.22495842D  02, 

2  0.74411568D  02,-0.414315760  03, -0 . 78754339D  02,  0.11254744D  02, 
2  0 .16021761D  03,-0.238621950  03, -0 . 50094687D  03, -0 . 68487854D  02, 
2  0.12254778D  02,-0.101619760  02, -0 .47219591D  01,  0.79729681D  01, 
2-0.21069574D  02,  0.22046490D  01,  0.B972B244D  01/ 

DATA  E/  0 .211031070  02, 

2-0 . 37959787D  03, -0 . 97489220D  02,  0.12900672D  03,  0.17949226D  02, 
2-0 . 12910931D  03,-0.557055740  03,  0.13524801D  02,  0.14696721D  03, 
2  0 . 17949528D  02,-0.329810140  00,  0.31028836D  02,  0.81657657D  01, 
2  0.22236961D  02,  0.39124892D  02,  0.81636799D  01/ 

Z-Vl 

DO  100  CTIM-1,2 

X- DREAL(Z) 

Y-DIMAG(Z) 

El 5-DCMPLX ( . 0D0 , .0D0) 

AB-CDABS (Z) 

IF (AB . EQ . 0 . ) GO  TO  90 

IF (X . GE . 0 .  .AND.  AB.GT.10.JGO  TO  80 

YA-DABS (Y) 

IF (X . LE . 0 .  .AND.  YA.GT.10.)GO  TO  80 

IF (YA-X.GE.  17.5  .OR.  YA.GE.  6 .5  .OR.X+YA.GE.5 .5  .OR.X.GE .  3 .  )GO  TO  20 
IF (X. LE . -9 . ) GO  TO  40 
IF (YA-X.GE. 2. 5)GO  TO  50 
IF (X+ YA.GE. 1.5) GO  TO  30 
10  N-6.+3.*AB 

E15-1 . / (N-l . ) -Z/N»*2 
15  N-N-l 

E15-1 . / (N-l . ) -Z*E15/N 
IF  (N .  GE .  3)  GO  TO  15 

E15»Z*E15-DCMPLX  ( .  577216+DLOG  (AB)  , DATAN2  (T,X)  ) 

GO  TO  90 
20  Jl-1 
J2-6 
GO  TO  31 

30  vfl**7 
J2-21 

31  8-DCHPLX(.0D0, .0D0) 

YS-Y*Y 

DO  32  I-Jl, J2 

XI- V(I)+X 
CF-W(I) / (XI*XI+YS) 

32  S-S+DCMPLX (XI*CF, -TA*CF) 

GO  TO  54 

40  T3-X*X-Y«Y 
T4-2 . *X*YA 
T5-X*T3-YA*T4 
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T6-X*T4+YA*T3  EXPJ.2 

OC-DCMPLX  (D  (11)  +D (12)  *X+D  (13)  *T3+T5-E  (12)  *YA-E  (13)  *T4, 

2  E (11) +E (12) *X+E (13) *T3+T6+D (12) *YA+D (13) *T4) 

VC-DCMPLX (D (14) +D  (15) *X+D (16) *T3+T5-E  (15) *YA-E (16) *T4, 

2  E (14) +E  (15) «X+E (16) *T3+T6+D (15) *YA+D (16) *T4) 

GO  TO  52 

50  T3-X*X-Y*Y 
T4-2 . *X*YA 
T5-X*T3-YA*T4 
T6-X*T4+YA*T3 
T7-X*T5-YA*T6 
T8-X*T6+YA*T5 
T9-X*T7-YA*T8 
T10-X»T8+YA*T7 

OC-DCMPLX  (D  (1)  +D  (2)  *X+D  (3)  *T3+D  (4)  »T5+D  (5)  *T7+T9-  (E  (2)  *YA+E  (3)  *T4 
2+B  (4)  «T6+E  (5)  *T8)  ,  B  (1)  +B  (2)  *X+E  (3)  *T3+E  (4)  *T5+B  (5)  *T7+T10+ 

3  (D  (2)  «YA+D  (3)  *T4+D  (4)  *T6+D  (5)  *T8) ) 

VC-DCMPLX  (D  (6)  +D  (7)  *X+D  (8)  *T3+D  (9)  *T5+D  (10)  *T7+T9-  (E  (7)  *YA+E  (8)  *T4 
2+E  (9)  »T6+E (10) *T8) ,E (6) +E (7) *X+E (8) *T3+E (9) *T5+E (10) *T7+T10+ 

3  (D  (7)  *YA+D  (8)  *T4+D  (9)  *T6+D  (10)  *T8)  ) 

52  EC-OC/VC 

S— EC/DCMPLX (X,  YA) 

54  EX-DEXP (-X) 

T— EX*DCMPLX  (DCOS  ( YA)  ,  -DSIN  (YA) ) 

E15-S*T 

56  IF  ( Y . LT .  0 .  )  E15-DCON JG  (E15) 

GO  TO  90 

80  E15— .409319/  (Z+ .193044) +.421831/  (Z+l .  02666)  +  .147126/  (24-2.56788)  + 

2 . 206335E-1/ (Z+4 . 90035) +.107401E-2/ (Z+8 . 18215) +.158 654E-4/ (2+ 

312 . 7342) + . 317031E-7/ (1+19 .3957) 

E15-E15*CDEXP(-Z) 

90  IF ( JIM.EQ . 1)W12— E15 
100  Z-V2 

Z-V2/V1 

TB-DATAN2  (DIMAG  <Z)  f DREAL  (Z)  )  -DATAH2  (DIMAG  (V2) , DREAL  (V2)  ) 

2+DATAN2  (DIMAG  (VI )  ,  DREAL  (VI ) ) 

AB-DABS (TB) 

IF (AB.LT.l.)TH-.O 
IF (TH.GT.l.)TH-6. 2831853 
IF  (TH.LT.-l.)TH— 6.2831853 
W12-W12-E15+DCMPLX ( . 0D0 ,  TH) 

RETORN 

END 
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SUBROUTINE  GRILL  (AX,  BAR,  DKD,  DKW, NEQ, NSD, NSV»,  TK,  VJ) 
C  GRILL  CALCULATES  THE  VOLTAGE  COLUMN  VJ(I)  rOR 
C  MONOPOLE  ON  CIRCULAR  DISK  IN  FREE  SPACE. 

IMPLICIT  REAL*8  (A-H) ,  (P-Z) 

COMPLEX* 16  EGZ,GM,GP,GI (20) ,VJ (1) , GII , QST, NST, VJ1 
DATA  PI, TP/3. 14159265359, 6. 28318530718/ 

IDM-20 

DO  20  I-1,NEQ 

20  VJ (I) -DCMPLX (.ODO, .0D0) 

VJ(1)-DCMPLX(1.D0, .0D0) 
ir (BAR . LE . 1 .) RETURN 
VJ (1)  -DCMPLX  ( .  0D0,  .0D0) 

DK-DKW 

SDK-DSIN(DK) 

CDK-DCOS (DK) 

BAL-DLOG (BAR) 

QST-DCMPLX ( . 0D0, 1 . / (4 . *BAL*SDK) ) 

BK-AK*BAR 

AKS«AK*AK 

BKS«BK*BK 

LIM-NSW+1 

IF (LIM . GT . IDM) LIM-IDM 
NPH**€ 

NPH-2* (NPH/2) 

NPP-NPH+1 
PHA-. 0174533*20. 

DPH-PHA/NPH 

PH-.O 

DO  90  LPH-1,2 
WST-DPH*QST/ (3.*PI) 

SGN—  1 . 

DO  80  IPH-l.NPP 
NF-3 . +SGN 
IF(IPH.EQ.1)NF-1. 

IF (IPH.EQ .NPP) WF"1 . 

CPH-DCOS (PH) 

IF ( IPH . GT . 1 ) GO  TO  40 
IF (LPH. GT. 1) GO  TO  40 
CPH-DCOS (DPH/10 . ) 

40  RS1-2.*AKS* (l.-CPH) 

RS2-AKS+BKS-2 . *AX*BK*CPH 
RH1-DSQRT  (RSI) 

RH2-DSQRT (RS2) 

CALL  CISI  (CA, CIN, SA, RH1 ) 

CALL  CISI2 (CB, CIN, SB,RH2) 

GI (1) -2 . *DCMPUC  (CB-CA, SA-SB) 

DO  50  1-2, LIM 
DZ-DK* (1-1) 

DZS— DZ*DZ 
RA-DSQRT (RS1+DZS) 

RB-DSQRT (RS2+DZS) 

CALL  CISI  (Cl, CIN, S1,RA+DZ) 

CALL  CISI  (C2,CIN, S2,RB+DZ) 

GP-DCMPLX (C2-C1, S1-S2) 

RAM-RS 1 / (RA+DZ ) 

RBM-RS2/ (RB+DZ) 

CALL  CISI  (Cl, CIN, SI, RAM) 

CALL  CISI  (C2 , CIN , S2 , RBM) 

GM— DCMPLX (C2-C1, S1-S2) 

EGZ-DCMPLX (DCOS (DZ) ,DSZN (DZ) ) 

50  GI (I) — GP*EGZ+GM/KGZ 

VJ(1)-VJ(1)+W*NST*  (GI  (2)  -CDK*GI  (1)  ) 

IF (NSW . LE . 1 ) GO  TO  78 
Kl-0 
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IA-NSD+1 

DO  60  I-IA,  NEQ 

Kl-Kl+1 

X2-X1+1 

X3-X2+1 

ir (K3 . GT . IDM) GO  TO  60 

GP-GI (Kl) -2 . *CDX*GI (X2) +GI (X3) 

VJ  < I ) -VJ ( I ) +WF*WST*GP 

CONTINUE 

SGN—SGN 

PH-PH+DPH 

DPH- (PI-PHA) /NPH 

PH-PHA 

CALL  CISI  (CA,CIN,SA,AX) 

CALL  CISI  (CB,CIN,SB,BX) 

R2-AK+DXD 
SR2— DSIN (R2) 

CR2-DCOS (R2) 

SDKD-DSIN (DKD) 

VI 1-  (SR2*  (CB-CA)  -CR2*  (SB-SA) )  /  (2 .  *BAL*SDXD) 
VJ(1)«VH+VJ(1) 

IF (NSD . LE . 1 ) RETURN 

V22-  (DSIN  (AK)  *  (CB-CA)  -DCOS  (AX)  *  (SB-SA)  )  /  (2 .  *BAL*SDXD) 
VJ  (2 )  -DCMPLX  (V22 ,  .0D0) 

RETURN 

END 
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SOBROOTINE  QDD  (CDKD,  SDKD,  SI,  S3,  Tl,  13,  TK,  IWZ,  NPH,  Z22) 

C  QDD  CALCULATES  Z22  -  MUTUAL  IMPEDANCE  OF  TMO  DISK  MODES. 
IMPLICIT  RZAL*8  (A-B) ,  (P-Z) 

COMPLEX*16  Z12, Z22, ZD 

DATA  PI, P3/3 . 14159265359, 9. 42477796077/ 

2  FORMAT (IX, 8F10 .2) 

5  FORMAT (1H0) 

6  FORMAT  (5X,  'PRINTODT  FROM  QDD') 

7  FORMAT <5X, 'DISK  DIPOLE  TO  DISK  DIPOLE') 

IF (IWZ . LE . 0) SO  TO  10 

WRITE  (17,  6) 

WRITE (17,7) 

WRITE  (17,5) 

10  PHA-. 0174533*2. 

PHB-. 0174533*20. 

NPH-2* (NPH/ 2) 

DPH-PHA/NPH 

PH-.O 

NPKP-NPH+1 
Z22-(.0D0, .0D0) 

DO  80  IPH-1,3 
ZD«( .0D0, .0D0) 

SGI— 1. 

DO  70  I-1,NPHP 

WF-3.+SGI 

IF(I.EQ.1)WF-1. 

IF(I.EQ.NPHP)WF-1. 

CPH-DCOS (PH) 

IF(I.EQ.l  .AND.  IPH.GT.DGO  TO  60 
C  NEXT:  DISK  DIPOLE  TO  DISX  DIPOLE. 

CALL  SKEWS  (SI, S3, Tl, T3,TK, CDKD, SDKD, CDKD, SDKD, CPH, Z12) 
60  PHD-57. 2 957 8 *PH 

ir  (IWZ .  GT .  0)  WRITE  (17,2) PHD, Z12 
SGI— SGI 
PH— PH+DPH 
70  ZD-ZD+WF*Z12 

Z22— Z22+DPH*ZD/P3 
PH-PKA 

DPH— (PHB-PHA) /NPH 
IF ( IPH . EQ . 1 ) GO  TO  80 
PH-PKB 

DPH-  (PI-PHB)  /NPH 
80  CONTINUE 

IF  (IWZ  .GT .  0)  WRITE  (17,5) 

RETORN 

END 
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SUBROUTINE  QDM  (AX,  DKD,  DKW,CDKD,  SDKD,  SDK,  31,  S3,  TK,  IWZ,NPH,  812) 
C  QDM  CALCULATES  Z12  -  MUTUAL  IMPEDANCE  OF  DISK  DIPOLE  AND  MODE 
IMPLICIT  R£AL*8  (A-H) ,  (P-Z) 

COMPLEX* 16  ZD, Z12, Z21, ZDM, ZM, P12, ZDD , PDM 
DATA  PI, P3/3 . 14159265359, 9 . 42477796077/ 

2  FORMAT <1X,8F10. 2) 

5  FORMAT (1H0) 

6  FORMAT <5X,' PRINTOUT  FROM  QDM') 

7  FORMAT (5X, 'DISK  DIPOLE  TO  DISK  MONOPOLE') 

6  FORMAT <5X, 'DISK  DIPOLE  TO  MIRE  MONOPOLE') 

IF (IMZ .LE . 0) GO  TO  10 
WRITE (17, 6) 

WRITE (17 ,7) 

WRITE (17,5) 

10  PHA-. 0174533*2. 

PHB-. 0174533*20. 

NPH-2* (NPH/2) 

HPHP-NPH+1 

DPH-PHA/NPH 

IDM-1 

IF (SI . CT . 10 . *AK) IDM-0 
PH-.O 

ZDD-(.0D0, .0D0) 

T2-AK+DKD 
DO  40  IPH-1,3 
ZD— ( . 0D0, .ODO) 

SGI— 1 . 

DO  30  I— 1 , NPHP 

WF-3.+SGI 

IF (I .EQ.l) WF-1 . 

IF<I.EQ.NPHP)WF-1. 

CPH-DCOS (PH) 

IF  (I . EQ .  1  .AND.  IPH.GT.DGO  TO  20 
C  NEXT:  DISK-DIPOLE  TO  DISX-MONOPOLE . 

CALL  ZSDM  (SI , S3, AX, T2, TK, CDKD, SDKD, SDKD, CPH, P12) 

20  PHD-57. 29578*PH 

IF (IWZ . GT . 0 ) WRITE ( 17 , 2) PHD,  P12 
SGI— SGI 
PH— PH+DPH 
30  ZD-ZD+WF«P12 

ZDD— ZDD+DPH*ZD/P3 
PH-PHA 

DPH- (PHB-PHA) /NPH 
IF ( IPH . EQ . 1 ) GO  TO  40 
PH-P8B 

DPH- (PI-PHB) /NPH 
40  CONTINUE 

IF  (IWZ  .GT .  0)  WRITE  (17, 5) 

W2— TK+DKW 

IF (IDM.EQ.O)GO  TO  100 
IF  (IWZ  .GT .  0)  WRITE  (17,  S) 

LPH-6 

LPH-2* (LPH/2) 

LPP-LPH-f  1 
PHA-. 0174533*20. 

DPH— PHA/LPB 
PH-.O 

ZDM— ( . ODO , .ODO) 

RMN-AK/100. 

DO  90  IPB-1,2 
PDM- (.0D0, .ODO) 

SGI— 1. 

DO  SO  I-1,LPP 


QDM.l 
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WF-3.+SGI 
Zr(Z.BQ.l)MF-l. 
IF(I.EQ.LPP)WF-1. 
CPH-DCOS  (PH) 
RH-AK*DSIN (PH) 
IF(I.EQ.1)RH-RMN 
AC-AK*CPH 
VI -SI -AC 


V3-S3-AC 

IF (I .EQ. 1  .AND.  IPH.GT . 1) GO  10  70 
NEXT:  DISX -DIPOLE  TO  WIRE-MONOPOLE. 

CALL  ZSDM  (V1,V3,TK,W2,RH,CDKD,SDXD,3DK, .0D0,F12) 
70  PHD-57. 29578*PH 

IF (IWZ.GT. 0) WRITE (17, 2) PHD,P12 
SGI— SGI 
PH-PH+DPH 
SO  PDM-PDM+WF*P12 

ZDM-ZDM+DPH*PDM/P3 
DPH- (PI-PHA) /LPH 
90  PH-PHA 

Z12— ZDM-ZDD 

IF (IWZ.GT. 0) WRITE (17, 5) 

RETURN 

NEXT:  DISK-DIPOLE  TO  WIRE -MONOPOLE.  _ 

100  CALL  ZSDM  (S1,S3,TK,W2,AX,CDKD,SDKD,SDK, .0D0,ZDM) 
Z12— ZDM-ZDD 
RETURN 
END 


I 

I 

I 

I 
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SUBROUTINE  QMM (AK, DKD, DKW, CDKD, SDKD, CDK, SDK, TK, IWZ,NPH, Zll) 
C  QMM  CALCULATES  Zll  -  SELF  IMPEDANCE  OF  MODE  1, 

C  WHICH  HAS  TERMINALS  AT  BASE  OF  KONOPOLE. 

IMPLICIT  REAL*8  (A-H) ,  (P-Z) 

COMPLEX* 16  FDM, FMD, FMM, PDH, PMD, PMM 
COMPLEX*  16  ZD,  ZDD,  ZDM, ZMD,  ZMM,  Z11,P11 
PATA  PI, P3/3 . 14 159265359, 9 . 42477796077/ 

1  FORMAT ( IX, 8F8.0) 

2  FORMAT (IX, 8F10.2) 

5  FORMAT (1H0) 

6  FORMAT (5X, 'PRINTOUT  FROM  QMM') 

7  FORMAT  (5X, 'DISK  MONOPOLE  TO  DISK  MONOPOLE') 

IF (IWZ .LE. 0) GO  TO  10 

WRITE (17 ,6) 

WRITE  (17, 5) 

10  AKS-AK«AK 
DKKP-TK+DKW 
CDKMP-DCOS (DKMP) 

SDKMP-DSIN (DKMP) 

ZMM-(.0D0, .0D0) 

ZMD-(.0D0, .0D0) 

ZDM- ( . 0D0, .0D0) 

LPH-6 

LPH-2* (LPH/2) 

LPP-LPH+1 

PHA- . 0174533*20 . 

DPH— PHA/LPH 
RMN-AK*DPH/10 . 

PH-.O 

DO  44  IPH-1,2 
FDM— ( . 0D0, .0D0) 

FMD— ( . 0D0 , .0D0) 

FMM— ( . 0D0, .0D0) 

SGI— 1. 

DO  40  I-l.LPP 

WF-3.+SGI 

IF ( I . EQ . 1 ) WF-1 . 

IF (I .EQ.LPP) WF— 1 . 

IF(I.EQ.l  .AND.  IPH.GT.DGO  TO  38 
CPH-DCOS (PH) 

8PH-DSIN (PH) 

DRG— 2 . * (1 . -CPH) 

R— AK*DSQRT (DRG) 

IF (I .EQ . 1) R-R+RMN 

C  NEXT:  WIRE  MONOPOLE  TO  WIRE  MONOPOLE. 

CALL  ZSMM  <-TK,DKW,  .  0D0 , DKW, R, CDKMP ,  SDKMP ,  SDK,  1 . DO ,  PMM) 
R“AK*SPH 

IF ( I. EQ.l) R-R+RMN 
Tl^AK* (1 . -CPH) 

T2-T1+DKD 

C  NEXT:  WIRE  MONOPOLE  TO  DISK  MONOPOLE. 

CALL  ZSMM  (-TK, DKW, T1,T2,R, CDKMP, SDKMP, SDKD,  .0D0,PMD) 

81- ^AK*  (1 . -CPH) 

82—  Sl+DKD 

C  NEXT:  DISK  MONOPOLE  TO  WIRE  MONOPOLE. 

CALL  ZSMM  (S1,82,TK,TX+DKW,R,CDXD,8DKD,SDK, . 0D0.PDM) 

38  FMM— FMM+WF  *PMM 
FMD— FMD+WF*PMD 
FDM— FDM+WF*PDM 
PHD— 57. 29578 *PH 

IF  (IWZ  .CT .  0)  WRITE  (17, 1) PHD, FMI,PMD,PDM 
SGI— SGI 
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40  PH-PH+DPH 

ZMM-ZMM+DPH*FMM/P3 

ZMD-ZMD+DPH*FMD/P3 

ZDM«ZDM+DPH*FDM/P3 

RMN-.O 

DPH-(PI-PHA) /LPH 
44  PH-PHA 

IF (IWZ. GT.O) WRITE (11,5) 

IF  (IWZ . GT .  0)  WRITE  (17,  7) 

PHA-. 0174533*2. 

PHB-. 0174533*20. 

NPB-2* (NPH/2) 

KPHP-NPH+1 

DPH-PHA/NPH 

PH-.O 

ZDD— ( . ODO,  .ODO) 

S2-AK+DKD 
DO  CO  IPH-1,3 
ED- (.ODO, .ODO) 

SGI— 1 . 

DO  50  I -1, HP HP 

WF-3.+SGI 

IF (I -EQ. 1) WF-1 . 

IF (I .EQ.NPHP) WF— 1 . 

CPB— DCOS (PH) 

IF (I .EQ. 1  .AND.  IPH.GT. l)GO  TO  48 
C  NEXT:  DISK  MONOP.OLE  TO  DISK  MONOPOLE. 

CALL  ESMM  <AK,  S2,  AK,  S2,  TK.CDKD,  SDKD,SDXD,CPB,P11) 
48  PHD— 57 .29578*PH 

IF  (IWZ .  GT .  0)  WRITE  (17,  2)  PHD,P11 
SGI— SGI 
PH-PH+DPH 
50  ZD— ZD+WF*P11 

ZDD-ZDD+DPH*ZD/P3 
PH— PHA 

DPH—  (PHB-PHA)  /NFH 
IF (IPH.EQ . 1) GO  TO  60 
PH-PBB 

DPH- (PI-PHB) /NPH 
60  CONTINUE 

Zll-ZDD-ZDM-ZMD+ZHM 
IF (IWZ. GT. 0) WRITE (17, 5) 

RETURN 

END 

C 

C 
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SUBROUTINE  SKEW (SI, S3, T1 , T3 , RHX, CDK, SDK, CDKT, SDKT, CPSZ, Z12) 
C  SKEW  CALCULATES  Z12  -  MUTUAL  IMPEDANCE  OF  CENTER-FED 

C  NONPLANAR-SKEW  SINUSOIDAL  DIPOLES  WITH  UNEQUAL  LENGTH. 

IMPLICIT  REALMS  (A-H) ,  (P-Z) 

COMPLEX* 16  Z12,EIN,EGDZ,CQX,EJXX,EM,EP 

COMPLEX* 16  E (2, 2) , F  (2, 2) , ESI , ES2 , ET1 , ET2, EXPA, EXPB, SGZI 
COMPLEX* 1 6  P11,P12,P21,P22 
DIMENSION  S  (3) ,  T  (3) 

DATA  ETA, PI, TP/376. 730366239, 3. 14159265359, 6.28318530718/ 

S(1)-S1 

S2-(Sl+S3)/2. 

S(2)-S2 
S  (3)  “S3 
T(1)-T1 
T2-(T1+T3) /2 . 

T(2)-T2 

T(3)«T3 

Z12“ ( . 0D0, .0D0) 

DPSI-CPSI 
CPSS-DPSI*DPSI 
IF (CPSS.GT.l.DO)CPSS-l.DO 
SPSI-DSQRT (1 .CO-CPSS) 

IF (DABS (CPSI) .LT. .999999)60  TO  10 

RHS-RHK*RHK 

RH2-SPSI* (T1+T3) / 2 . 

RHS-RHS+RH2*RH2 

SGN-1. 

IF (CPSI .GT. .0) GO  TO  BO 
SGN—  1. 

T(l)— T3 
T(2)— T2 
T(3)— T1 
GO  TO  80 
10  D-RHK 
DSQ-D«D 
CD-D/SPSI 
BD-CD*DPSI 
EB-DEXP(-BD) 

EC-DEXP (-CD) 

CST— ETA/ (16 . *PI*SDK«SDKT) 

TA-T1 

TB-T2 

DO  70  ITT-1,2 

IF  ( ITT .  EQ .  1 )  ET1— DCMPLX  (DCOS  (TA)  ,DSIN(TA) ) 
IF(ITT.EQ.2)BT2-DCMPLX(DCOS(TB) ,DSXN(TB)) 

TD1-TA 

TD2-TB 

TS1-TD1*TD1 

TS2— TD2*TD2 

SA-S1 

SB-82 

DO  60  ISS-1,2 

IF  (ISS.EQ. 1)KS1— DCMPLX (DCOS (SA) ,DSZN(SA) ) 

IF  (ISS . EQ . 2 ) ES2-DCMPLX (DCOS (SB) , DSXM (SB) ) 

DO  20  X-1,2 
DO  20  L— 1, 2 
20  E(K,L)-(.0D0, .0D0) 

SI— SA 

DO  50  1-1,2 

FI—  (-1) **I 

SDI-SI 

SIS— SDI*8DI 

ST1-2 . *SDI*TD1*DFSI 
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ST2-2 . *SDI*TD2*DPSI 
Rl-DSQRT (DSQ+SIS+TS1-ST1) 

R2-DSQRT (DSQ+SIS+TS2-ST2) 

EK-EB 

DO  40  K-1,2 
FK-(-l)**K 
SK-FK*SDI 
BL-EC 

DO  30  L-1,2 

FL-(-l)**L 

EKL-EK*EL 

XX-FK*SD+FL*CD 

TLl-rL*TDl 

TL2-FL*TD2 

RR1-R1+SK+TL1 

RR2-R2+SK+TL2 

CKLL  EXP  J  (DCMPLX  (XX,RR1)  ,  DCMPLX  (XX,  RR2)  ,EXPA) 

CALL  EXPJ  (DCMPLX ( -XX, RR1)  , DCMPLX ( -XX, RR2)  ,EXPB) 

E  (K,  L)  -E  (K,L)  +FI*  (SXPA*EKL+EXPB/EKL) 

30  EL-1 . /EC 
40  EK-l./EB 

IF (I . EQ . ISS) GO  TO  50 

ZD— SDI*DPSI 

ZC-ZD 

EGZ I -DCMPLX (DCOS (ZC) ,DSZN(ZC) ) 

RR1— R1+ZD-TD1 
RR2-R2+ZD-TD2 

CALL  EXPO  (DCMPLX  (.0D0,RR1),  DCMPLX  (.0D0,RR2>  ,EXPB) 

RJU-Rl-ZD+TDl 

RR2-R2-ZD+TD2 

CALL  EXPJ  (DCMPLX <.0D0,RR1)  .DCMPLX (.OD0.RR2)  ,EXPA) 

F  (1, 1>-  (  .  CD0, 2  .DO)  *SDK*EXPA/EGZI 
T (1, 2) - ( . 0D0, 2 .DO) *SDX*EXPB*EGZI 
50  SI— SB 

IF (ITT. EQ. 2) GO  TO  54 

IK  (ISS  .EQ . 1) P22-CST* ( (F (2, 1) +B (2, 2) *ES1-E (1,2) /ESI) *ET1 
A+  (-F  (2, 2)  -E  (2, 1)  *ES1*E  (1,1)  /ESI)  /ET1) 

IF  (ISS .EQ . 2) P12-CST* ( (-F (1, 1) -E (2 , 2) *ES2+E (1,2) /ES2) «ET1 
B+(  F  (1, 2)  +E  (2,1)  -ES2-E  (1,1)  /E52)  /*T  1 ) 

GO  TO  58 

54  IF (ISS .EQ . 1) P21-CST*  (  (-F  (2, 1) -E (2, 2) *ES1-HB (1, 2) /ESI) *BT2 
C+(  F (2, 2) +E (2, 1) *ES1-E (1, 1) /ESI) /ET2) 

IF  (ISS.EQ.2)P11— CST* ( (r (1 , 1) +E (2, 2) *ES2-E (1, 2) /ES2) *ET2 
D+  (-F  (1, 2)  -E  (2, 1)  *ES2+E  (1, 1)  /ES2)  /ET2) 

58  SA-S2 
SO  SB— S3 
TA-T2 
70  TB-T3 

Z12-P11+P12+P21+P22 

RETURN 

80  DO  100  1-1,3 
SI-S(I) 

CI-1. 

IF  (I  .EQ.2)  Cl— 2  .  *CDK 
CQX— ( . 0D0 , .0D0) 

DO  90  J-1,3 
TJ-T ( J) 

CJ-1. 

IF ( J.EQ . 2)  CJ— 2 . *CDKT 
DE-TJ-SI 

R-DSQRT (RHS+DZ«DZ) 

X-R+DZ 

IF  (DZ . LT .  .  0)  X-RHS/  (R-DZ) 

CALL  CISI  (COS I , CIM , 8 INI ,  X) 

EP— DCMPLX (COS I, -SIM!) 

X-R-DZ 
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IF  (DZ  .  ST .  .  0)X-RBS/  (R+DZ) 

CALL  CISI  (COSI, CIN, SINI,X) 
EM-DCMPLX (COSI, -SINI) 

EGDZ-DCMPLX  (DCOS  (DZ) ,  DSIN  (DZ)  ) 

90  CQX-CQX+CJ* <XP*EGDZ+EM/EGDZ) 

X00  Z12—Z12+CI*CQX 

Z12-  SGN*ETA»Z12/ (8 . *PI*SDK*SDKT) 

RETURN 

END 
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SUBROUTINE  SKEWS (SI, S3, T1 , T3, RHK,CDK, SDK, CDKT, SDKT, CPSI, Z12) 
C  SKEWS  CALCULATES  Z12  -  MUTUAL  IMPEDANCE  OF  CENTER-FED 
C  COPLANAR - SKEW  SINUSOIDAL  DIPOLES  WITH  UNEQUAL  LENGTHS. 

IMPLICIT  REAL* 6  (A-H) ,  (P-Z) 

COMPLEX  * 1 6  Z12,EIN,EGDZ,CQX,EJXX,EM,EP 
DIMENSION  S  (3)  ,  T  (3)  • 

DATA  ETA, PI, TP/376.730366239, 3.14159265359, 6.28318530718/ 
8(1) -SI 

S (2)  —  (S1+S3) /2 . 

S  (3)— S3 
T(1)*T1 

T  (2)  —  (T1+T3)  /2 . 

T  (3)-T3 

Z12— ( . 0D0, .0D0) 

DPSI-CPSI 

IF (DABS (CPSI) .LT. .999999) GO  TO  10 

RHS-RHK*RHX 

CPSS-DPSI*DPSI 

IF (CPSS .GT.1.D0) CPSS-1.D0 

SPSI-DSQRT (1 .DO-CPSS) 

RH2-SPSI* (T1+T3) /2 . 

RHS— RHS+RH2*RH2 
8GN-1. 

IF (CPS1 . GT . .0) GO  TO  60 
SGN — 1 . 

T(l)—  T3 
1  (2)  — T  (2) 

T  (3)  — T1 
GO  TO  60 
10  DO  50  1-1,3 
SI-S(I) 

SIS-SI*SI 

CI-1. 

IF  ( I .  EQ .  2 )  Cl— 2 .  *CDK 
DO  50  J-1,3 
TJ-T  ( J) 

TJS-TJ*TJ 

R-DSQRT (SIS+TJS-2 . *SI*TJ*DPSI) 

CJ-1. 

ir  ( J .  EQ .  2 )  CJ— 2  .  *CDKT 
CQX— ( . 0D0, . 0D0) 

DO  40  K-1,2 
FK-(-l.)**X 
DO  40  L-1,2 
FL— (-1 . ) **L 
XXD— FK*SI+FL*TJ 
XX-XXD 

E JXX-DCMPLX  (DCOS  (XX)  ,  DSIN  (XX)  ) 

XXX-R+XXD 
X-DABS (XXX) 

CALL  CISI  (COSX , CXM , SINI ,  X) 

IF  (XXX . LT .  .  0)  SINI— SINI 
CQX— CQX+DCMPLX  (COSI,  -SINI)  *EJXX*FX*FL 
40  CONTINUE 

Z12-Z12+CQX«CI*CJ 
50  CONTINUE 

Z12— ETA*Z12/ (6 . *PI*SDK*8DKT) 

RETURN 

60  DO  80  1-1,3 
81— S  (I) 

CI-1. 

IF  (I. EQ. 2)  Cl— 2.*CDK 
CQX-(.0D0, .0D0) 
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DO  70  J-l#3 
TJ-T  ( J) 

CJ-1. 

IF  ( J .  EQ  •  2 )  CJ—  2 .  *CDKT 
DZ-TJ-SI 

R-DSQRT (RHS+DZ*DZ) 

X-R+DZ 

ir  <DZ . LT . . 0) X-RHS/ (R-DZ) 

CALL  CISI  (COSI/ CIN, SINI, X) 
EP-DCMPLX (COSI , -SINI) 

X-R-DZ 

ir(DZ.GT. . 0 ) X“RHS / (R+DZ ) 

CALL  CISI  (COSI, CIN, SINI,  X) 
EM-DCMPLX (COSI, -SINI) 
EGDZ-DCMPLX  (DCOS  (DZ),DSIN  (DZ)  ) 

70  CQX-CQX+CJ* (EP*EGDZ+EH/EGDZ) 

80  Z12-Z12+CI*CQX 

Z12-SGN*ETA*Z12/ (8 . *PI*SDK*SDKT) 

RETURN 

END 
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SUBROUTINE  8 KENT (AX, SI, S3, Tl, T3, CDK, SDK, CDKD, SDKD, INZ, Z12) 
C  SKENT  CALCULATES  Z12  *  MUTUAL  IMPEDANCE  OF  MIRE  DIPOLE 
C  AND  DISK  DIPOLE. 

IMPLICIT  REAL*8  (A-B) ,  (P-Z) 

COMPLEX*]. 6  P12,  Z12,QX2 

DATA  PI, P3/3. 14159265359, 9. 42477796077/ 

2  FORMAT (IX, 8F10.2) 

5  FORMAT (1H0) 

6  FORMAT (5X, 'PRINTOUT  FROM  SKEWT') 

7  FORMAT (5X, 'MIRE  DIPOLE  TO  DISK  DIPOLE') 

IF (INZ .LE . 0) GO  TO  20 

MRITS  (17,  6) 

MRITE (17, 7) 

MRITE (17 , 5) 

20  RMN-AK/100 . 

NPHW6 

NPH-2* (NPH/2) 

NPP-NPH+1 

PHA-. 0174533*20. 

DPH-PHA/NPH 
Z12-(.0D0, .0D0) 

PH-.O 

C  MIRE  DIPOLE  TO  DISK  DIPOLE. 

DO  80  1PB-1,2 
Q12-(.0D0, .0D0) 

SGI— 1. 

DO  60  1-1, NPP 

NF-3.+SGI 

IF ( I . EQ . 1 ) NF— 1 . 

IF(I.EQ.NPP)NF-1. 

RH-AK*DSIN (PB) 

IF ( I . EQ . 1 ) RH-RH+RMN 
AC-AK*DCOS (PH) 

V1-T1-AC 

V3-T3-AC 

IFd.EQ.l  .AND.  IPH.EQ.2) GO  TO  50 

CALL  SKEN  (SI, S3, VI, V3,RB, CDK, SDK, CDKD, SDKD,  .0D0,F12) 

50  QI2-Q12+MF*P12 

PHD-57. 2 957 8 *PH 
IF  ( IMZ .  GT .  0 )  MRITE  ( 17 , 2 )  PHD ,  P12 
SGI— SGI 
60  FH-PH+DPH 

Z12-Z12+DPH*Q12/P3 

RMN-.O 

DPH— (PI-PBA) /NPH 
80  PH-PHA 

IF(IMZ.GT.0)MRITE(17,5) 

RETURN 

END 


SKEWT 
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SUBROUTINE  SPART (AK, DKD, DKW, MAX, IWZ, Z, Zl)  SPART.l 

C  SPART  CALCULATES  Z  -  MUTUAL  IMPEDANCE  OF  TWO  WIRE  DIPOLE  MODES, 

C  AND  Z1  -  MUTUAL  IMPEDANCE  BETWEEN  A  WIRE  DIPOLE  MODE  AND  MODE  1. 

IMPLICIT  REAL*  8  (A-H)  ,  (P-Z) 

COMPLEX*  1 6  EID  (20)  ,EM(20)  ,EP  (20) ,  Z  (I)  ,  ZI  (1) 

COMPLEX* 16  CEM, CEP, EMD, EPD, EMD2, EPD2, Zll, Z22 , Gil , Qll 
DIMENSION  CID  (20)  ,  SID  (20)  ,CM(20)  , CP  (20)  ,  SM (20)  ,  SP  (20) 

DATA  GAM, P2/. 577215664, 1.57079632/ 

DATA  ETA, PI/376. 727, 3. 14159/ 

IDM-20 

1  FORMAT (3X, 'MUST  INCREASE  DIMENSIONS  IN  SUBROUTINE  S PART') 

2  FORMAT (3X, 'ACTUAL  DIMENSION  IDM  -  ',I5,6X, 

2 'REQUIRED  DIMENSION  MAX2  -  ',15) 

3  FORMAT(lX, 8ri0.2) 

5  FORMAT (1H0) 

6  FORMAT (5X, 'PRINTOUT  FROM  SPART') 

7  FORMAT  (5X, 'FIRST:  WIRE  DIPOLE  TO  MODE  ONE') 

8  FORMAT (5X, 'THEN:  WIRE  DIPOLE  TO  WIRE  DIPOLE') 

IF (MAX. LE.0) RETURN 

MAX2-MAX+2 
DO  14  1-1 , MAX 
Zl  (!)•(. 0D0,  .0D0) 

14  Z(I)-(.0D0, .0D0) 

IF (MAX2 . LE . XDM) GO  TO  16 
WRITE (17, 1) 

WRITE (17, 2) IDM,MAX2 
RETURN 
16  DK-DKW 

IF (IWZ .LE . 0) GO  TO  18 
WRITE  (17, 6) 

WRITE (17,7) 

WRITE  (17,  8) 

WRITE(17,5) 

18  TDK-2. *DK 

SDKD-DSIN (DIO) 

S11-.0 

SI 3— TDK 

S21-DK 

S23— 3 . *DK 

DO  20  N-1,MAX2 

I-N-l 

DZ— I*DK 

CID  (N)-DCOS  (DZ) 

SID  (N)-DSIN  (DZ) 

20  KID  (N)  — DCMPLX  (CID  (N)  ,  SID  (N) ) 

CDK-DCOS (DK) 

SDK— DSIN (DK) 

EPD— DCMPLX (CDK, SDK) 

KMD-DCMPLX (CDK, -SDK) 

KPD2— EPD*EPD 
ZMD2-EMD*EKD 
CBM— 2 . *CDK+EMD 
CEP-2 . *CDK+EPD 
AX2— AK*AK 

CSS-ETA/ (8 . *PI*SDK*SDK) 

||p  Ha  5 

NPH-2* (NPH/2) 

NPP-NPH+1 
PHA-. 0174533*20. 

DPB— PHA/NPB 
PH-.O 

DO  100  JPH-1,2 

CST— DPH*ETA/  (24 . *PI*PI*SDK*SDK) 
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C22-DPB/(3.4pi) 

SGI,— 1 . 

DO  80  IPB-1,NPF 
CPH-DCOS (FH) 

SFB-DSIN (FB) 

IF (IPH.GT.l)GO  TO  30 
IF (JFB.GT.l)GO  TO  30 
FB0-DPB/10. 

CPH-DCOS (FHO) 

SFH-DSZN (FBO) 

Tl^AK® (l.-CFH) 

T2-T14DKD 

RH-AK*SPH 

DRG— 2 . *AK2* (1 . -CPB) 

RK-DSQRT (DRG) 

rs-drg 

WF-3.4SGN 
IF  (IPH.EQ . 1)MF- 1 . 

IF (ZFB . EQ . NFP) KF-1 . 

BST-KF*CST 

M22«Mr»C22 

DO  40  N-1,MXX2 

I-N-l 

DZ— I*DK 

DZS— DZ*DZ 

R-DSQRT (RS+DZS) 

ARG-R4DZ 

IF(N.EQ.l)  ARG— RK 

CALL  CISI (CP (N) , CIM, SF (M) , ARG) 

EP  (H)  — DCMPLX  (CP  (N)  ,  -SP  (K)  ) 

IF (N . GT . 1 ) GO  TO  38 
CM  (1)  —CP  (1) 

SM(1)-SP(1) 

EM(1)-EP(1) 

GO  TO  40 
ARG— RS/ ARG 

CALL  CISI(CM(H),Cm,SH(H),ARG) 
EM  (N)  -DCMPLX  (CM  (N) ,  -SM  (N)  ) 
CONTINDE 


R-4  .  •  (-CM (2)  +2 .  *CP  (1)  -CP  (2)  ) 

A+2 .  *CID  (3)  *  (+CH (3)  -2 .  *CM  (2)  42 . *CP  (1)  -2 .  *CP  (2)  +CP  (3)  ) 
B+2  .  *SID  (3)  *  (-SM (3)  +2 . *SM  (2)  -2 .  *$P  (2)  4SP  (3)  ) 

X— 4 .  *  (SM(2)  -2 .  *SP  (1)  +SP  (2) ) 

C+2 .  *CID  (3)  *  (-SM  (3)  42 .  *5M  (2)  -2 .  *SP  (1)  42 .  *SP  (2)  -SP  (3) ) 
D+2 .  *SID  (3)  *  (-CM  (3)  +2 .  *CM  (2)  -2 .  *CP  (2)  4CP  (3) ) 

Q1 1-CSS ‘DCMPLX (R,X) 

Z  ( 1 )  — Z  (1 )  4W3T*DCMPLX  (R,  X) 

ZU— 6 .  *CDK*KP  (1)  42 .  *  (EPD4CDK)  **P  (2)  -EPD*EP  (3)  4 
22 .  *  (EMD4CDX)  *KM (2)  -EMD*EM (3) 

CALL  ZSDM(S11,S13,T1,T2,RB,CDX, SDK, SDKD, .0D0,Z22) 
Gil— CSS*Z11-Z22 


El (1)-Z1 (1) 4WST*Z11-W22«Z22 
PBD— 57 . 29578*PB 


IF(IWZ.GT.O)WRITE(17,3)PBDrGll,011 
IF (MAX . EQ . 1 ) GO  TO  70 

R-2 . *CID (2)  *  (*CM(3)  43.*CM(2)-4.*CP (1)43. *CP (2) -CP  (3) ) 
E42 .  *SID  (2)  *  (48M (3)  -2 .  *SM  (2)  42 . «SP  (2)  -8P  (3) ) 

F4CID (4)*  (4CM (4) -2 . *CM (3) 4CM (2) 4CP (2) -2 . *CP (3) 4 CP (4) ) 
G4SID  (4)  •  (-8M  (4)42.  *8M  (3)  -8M  (2)  4SP  (2)  -2 .  *8P  (3)  48P  (4)  ) 
X-2 .  *CID  (2 )  •  (EM  (3)  -3 .  *SM  (2)  44 .  *SP  (1 )  -3 .  «8P  (2)  4  8P  (3)  ) 
B  42 .  *SID  (2)  *  (CM  (3)  -2 .  *CM  (2)  42 .  *CP  (2)  -CP  (3) ) 

I4CZD  (4)  *  (-8M  (4)  42 .  *8M(3)  -SM  (2)  -SP  (2)  42 .  *8P  (3)  -SP  (4)  ) 
J4SID  (4)*  (-CM  (4)42.  *CH  (3)  -CM  (2)  4CP  (2)  -2 .  «CP  (3)  4CP  (4)  ) 

Z  (2)  -Z  (2)  4WST*DCMFLX(R,X) 

Zll-2 .  *EP  (1)  4CZM*  (ZMD*SM  (3)  -EPD*EP  (2) ) 4 
2 CEP*  (KPD*EP  (3)  -EMD*EM (2)  )  -EPD2*EP  (4)  -EMD2*EM  (4) 
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ESDM(S21,S23,T1,T2,RH,CDK, SDK, SDKD,  .0D0,822) 
81  (2)  -81  (2)  +HST*Z11-H22*822 
IF (MAX. EQ. 2) GO  TO  70 
81— DK 

DO  80  N-3 , MXX 

Ml-N-1 

M2-N-2 


Nl-N+1 

i  AvTfN  (ini  4- EM  (Ml)  /BIO  (Ml)  +CEP*  (BP  (HI)  ‘BID  (H)  -EM  (H)  /BIO  (H)  ) 
wKS  -“(“)  /5£  (H) >  -HP  (H2)  ‘BIO  (HI)  -BM(H2)  /BIO  (HI) 

81— Sl+DlC 

ESDM{Sl,S3,Tl,T2,RH,CDK,SDKfSDKD,  .000,822) 

81 (H) -81 (H) +MST*811-W22*822 
CPA-CP  (M2)  -2 .  *CP  (Ml)  +CP  (H) 

CFB— 2 .  *CP  (H)  -CP  (Ml)  -CP  (HI) 

CPC-CP  (H2 )  -2 .  *CP  (HI )  +CP  (H) 

CMA-CM (M2)-2.*CM(M1)+CM(H) 

CMB— 2 .  *CM  (H)  -CM  (HI )  -CM  (Ml) 

CMC— CM (H2)  -2 .  *CM (HI)  +CM  (H) 

SPA-SP  (M2)  -2 .  *SP  (Ml)  +SP  (H) 

SPB— 2 .  *SP  (H)  -8P  (Ml)  -8P  (HI) 

8PC-SP  (H2)  -2 .  *3P  (HI)  +SP  (H) 

SMA-SM (M2)  -2 .  *SM  (Ml)  +SM (H) 

SMB-2 .  <*SM (H)  -SM (Ml)  -8M  (Ml ) 

SMC-SM  (N2)  -2 .  *SM  (HI)  +SM (H)  mmi 

R-CID  (M2)  *  (CPA+CMA)  *2 .  ‘CIO  (H)  *  (CPB+CMB)  +2 .  *SID  (H)  (SPB-SMB) 

K  +CID  (H2)  *  (CPC+CMC)  +SID  (H2)  *  (SPC-SMC) 

KifoJ!  ’4  .'"«??•  <»»„»>«.  *sn >.»>  *  «*.-«»> 

t  -CID  (H2)  *  (SPC+SMC)  +SID  (H2)  •  (CPC-CMC) 

IF  (H .GT . 3)  X-X+SIO  (M2)  *  (CPA-CMX) 

60  8  (H)  —8  (H)  +HST*DCMPLX  (R,  X) 

70  PH-PH+DPH 
80  SGN--SGN 

DPH— (PI-PHA) /HPH 
100  PH-PHA 

IF  (IHZ . GT .  0)  WRITE  (17, 5) 

RETORN 

END 


C 

c 


A-43 


ou 


39 

SUBROUTINE  ZSDM(S1,S3,T1,T2,RHK,CDK, SDK, SORT, CPSX, Z12) 

C  CALCULATES  Z12  -  MUTUAL  IMPEDANCE  BETWEEN  SINUSOIDAL  DIPOLE 

C  AND  SINUSOIDAL  MONOPOLE  KITH  SKEW  ORIENTATION. 

IMPLICIT  REAL*8  (A-H)  ,  (P-Z) 

COMPLEX* 16  E (2,2) ,F(2,2),ES2,ET2,EXPA,EXPB,EGZI.ES1 
COMPLEX*l6  CQX,  E JXX,  Z12 ,  EP1 ,  EMI , Pll , P21 , SS)ZtEM, EP 
DIMENSION  S (3) 

DATA  ETA, PI,  TP/376 . 730366239, 3 . 14159265359, 6.29318530718/ 

S(1)-S1 

S2-<Sl+S3)/2. 

S(2)-S2 
S (3) “S3 

Z12-(.0D0, .0D0) 

DPSX-CPSX 
CPSS-DPSX*DPSX 
IF (CPSS . GT . 1 . DO ) CPSS-1  .DO 
SPSI-DSQRT (1 .DO-CPSS) 

IF (DABS (CPSI) .LT.. 999999) GO  TO  10 

RHS-RHK*RHK 

RH2-3PSI* (T1+T2) /2 . 

RHS-RHS+RH2*RH2 

8GN-1. 

IF (CPSI. GT. .0)GO  TO  80 

S  (1)  ■'“S3 

S(2)--S2 

S(3)— SI 

SGN— 1 . 

GO  TO  80 
10  D-RHK 
DSQ-D*D 
CD-D/SPSI 
BD-CD*DPSI 
EB-DEXP (-BD) 

EC-DEXP (-CD) 

TD1-T1 

TD2-T2 

TS1-TD1*TD1 

TS2-TD2*TD2 

CST-- ETA/ (16 . *PI*SDK*SDKT) 

SA-S1 

SB-S2 

ET2-DCMPLX  (DCOS  (T2)  ,D3XN  (T2)  ) 

DO  60  ISS- 1,2 

IF  (ISS  .EQ. 1) ESI— DCMPLX (DCOS  (SA)  ,DSIN(SA)) 

IF  (ISS .  EQ .  2)  ES2-DCMPLX  (DCOS  (SB)  ,  DSIN  (SB)  ) 

DO  20  K-1,2 
DO  20  L-1,2 
20  K(K,L)-(.0D0, .0D0) 

8I-SA 

DO  SO  1-1,2 

PI-(-l)**I 

EDI— SI 

SIS— 5DI *5DZ 

8T1-2 . *SDI*TD1*DPSI 

8T2-2 . •SDI*TD2*DPSI 

Rl-DSQRT (DSQ+SIS+TS1-8T1) 

R2-DSQRT (DSQ+SIS+TS2-8T2) 

EX-EB 

DO  40  K-1,2 
FK-(-l) **K 
8K-FK*SDI 
EL-EC 

DO  30  L-1,2 


ZSDM.l 
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fl-(-d**l 

EXL-EK*BL 

XX-FK*BD+rL*CD 

TL1-FL*TD1 

TL2«FL*TD2 

RR1-R1+SK+TL1 

RR2-R2+SK+TL2 

CALL  EXPJ  (DCMPLX  (XX,  RJU)  , DCMPLX (XX, RR2)  , EXPA) 
CALL  EXPJ  (DCMPLX  (-XX,RR1) , DCMPLX ( -XX, RR2)  , EXPB) 

E  (K,  L)  -E  (K,  L)  +FI*  (EXPA*EXL+EXPB/EKL) 

30  EL-1 . /EC 
40  EX-l./EB 

XF (X.EQ.XSS)GO  TO  50 

ZD— SDI*DPSX 

ZOZO 

EGZI— DCMPLX  (DCOS  (ZC)  ,DSIN(ZC) ) 

RR1— R1+ZD-TD1 
RR2-R2+ZD-TD2 

CALL  EXPJ  (DCMPLX  ( .0D0,RR1) , DCMPLX ( .  0D0,RR2)  ,EXPB) 

RRI-R1-ZD+TD1 

RR2— R2 -ZD+TD2 

CALL  EXPJ  (DCMPLX  ( . ODO,RR1)  r  DCMPLX  ( . ODO,RR2)  ,EXPA) 
F  (1, 1)  -  (.  ODO,  2  .DO) *  SDK*EXPA/EGZX 
F (X, 2) - ( . ODO , 2 .DO) *SDK*EXPB*EGZI 
50  SI— SB 

XF(ISS.EQ.l) 

AP21-CST*  { (-F  (2, 1)  -E  (2, 2)  *ES1+B  (1,2)  /ESI)  *ET2 
B+  (  F  (2, 2)  +E  (2,  1)  *ES1-S  (1,1)  /ESI)  /ET2) 

IF  (ISS.EQ.2) 

CPU— CST*  (  (  F  (1, 1)  +E  (2,2)  -ES2-E  (1,2)  /BS2)  *ET2 
D+  (-F  (1, 2)  -E  (2, 1)  *ES2+E  (1,1)  /ES2)  /ET2) 

SA-S2 
60  SB— S3 

Z12-P11+P21 

RETURN 

80  DO  100  1-1,3 
CX-1. 

IF  (X  . EQ.2)  Cl— 2  .*CDK 
SI-S(X) 

TJ-T1 

DO  90  J-1,2 
DZ-TJ-SI 

R-DSQRT <RHS+DZ*DZ) 

X-R+DZ 

XF (DZ . LT . . 0 ) X-RHS/ (R-DZ) 

CALL  CXSX(COSX,CXN,SXia,X) 

EP— DCMPLX (COSX, -SXNX) 

X-R-DZ 

IF  (DZ  .CT . .  0)  X-RHS/  (R+DZ) 

CALL  CXSX (COSX, CXN, SXNX, X) 

EM-DCMPLX (COSX , -SXNX) 

XF ( J.EQ.2) GO  TO  90 

EP1-ZP 
EMI— EM 

90  TJ-T2 
X-T2-SI 

EGDZ-DCMPLX  (DCOS  (X)  ,DSXN  (X)  ) 

Z12— Z12+CI*  (  (EP-EP1)  »EGDZ+  (EM-EK1)  /ECDZ) 

100  CONTINUE  _ _ 

Z12-  SGN*ETA* Z12 /  (• . *PI*SDK*SDKT) 

RETURN 

END 


C  U 

c 

SOB  ROUTINE  ZSMM (SI, S2,T1, T2,D, CDS, SDS,  SDT,CPSI,P11)  ZSHM. 

C  CALCULATES  MUTUAL  IMPEDANCE  OF  COPLAKAR-SKEW 

C  SINUSOIDAL  MONOPOLES . 

IMPLICIT  REAL* 8  (A-H)  ,  (P-Z) 
complex* 16  E<2,2) ,r(2,2) ,gam,pii,pi2,p21,p22 
complex* 16  egzi,esi,es2,eti,et2,expa,expb 
complex* i 6  egdz,em,ep,emi,epi 

DATA  ETA, GAM, PI/376. 730366239, (.0D0,1.D0) ,3.14159263359/ 

DD-D 

DPQ— DD*DD 

DPSI-CPSI 

CPSS-DPSI*DPSI 

IF (CP SS . GT . 1 . DO ) CPSS-1 .DO 

8PSI-DSQRT(1.D0-CPSS) 

8GDS-SDS 

IF  (S2 .  LT .  SI)  SGDS— SDS 
SGDT-SDT 

IF (T2 .LT .Tl) SGDT-- SDT 

IF (DABS (CPSI) .LT . .999999)00  TO  € 

DO-SPSI* (T1+T2) /2 . 

DSQ“DPQ+DO*DO 
GO  TO  110 

6  ESl-DCMPLX(DCOS  (SI)  ,DSIN(S1) ) 

ES2-DCMPLX (DCOS (S2) ,DSIN(S2) ) 

ET1-DCMPLX (DCOS  (Tl)  ,DSIN  (Tl)  ) 

ET2-DCMPLX  (DCOS  (T2)  , OSIN  (T2)  } 

ID-1 

IF  (D .  EQ . .  0)  ID— 0 

TD1-T1 

TD2-T2 

CD-DD/SPSI 

C-CD 

BD— CD*DPSI 
B-BD 
EB-.O 
EC-.O 

ir (ID.EQ.O)GO  TO  8 
EB-DEXP (-B) 

EC-DEXP (-C) 

8  DO  10  K-1,2 
DO  10  L-1,2 
30  K (K, L) — ( . 0D0 , . 0D0) 

TS1— TD1*TD1 
TS2— TD2*TD2 
SI-SI 

DO  100  1-1,2 

FI-(-l)**I 

SOI— SI 

SIS— SDI*SDI 

8T1-2 . *SDI*TD1*DPSI 

ST2-2 . *SDI *TD2»DPSI 

Rl-OSQRT  (DPQ+SIS+TS1-ST1) 

R2-DSQRT (DPQ+SIS+TS2-ST2) 

EK-XB 

DO  50  K-1,2 
FK-(-l)**K 
SK-FK*SDI 
EL-EC 

DO  40  L-1,2 
FL-(-l)**L 
EKL— EK*EL 
XX— FK*BD+rL*CD 
TL1“FL*TD1 
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TL2«FL*TD2 
RRI-R1+SK+TL1 
RR2-R2+SK+TL2 
AXX-DABS (XX) 

ZF (AXX.GT.DABS (RRZ) /100 . ) GO  TO  28 
ZF (AXX. GT.DABS  (RR2) /100. )G0  TO  28 
ZF(AXX.GT. .001)60  TO  28 
ZF (RRZ/RR2 . LT . . 0) GO  TO  28 
CALL  CZSZ (COSZ, CZN,SZNZ,RR1) 

CALL  CZSZ (COS2, CZN, SZN2, RR2) 

EXPA-DCMPLX (COS2-COS1, SIH1-SIM2) 

E  (K,  L)  -E  (K,  L)  +FZ  *EXPA*  (EKL+1 .  /EEL) 

GO  TO  40 

28  CALL  EXP  J  (DCMFLX  (XX, RR1)  , DCMPLX  (XX,  RR2)  , EXPA) 
CALL  EXPJ (DCMPLX ( -XX, RRZ)  , DCMPLX (-XX, RR2)  ,EXPB) 
E  (K,  L)  -E  (K,  L)  +FI*  (EXPA*EXL+EXPB/EKL) 

40  EL-Z./EC 
50  EK-Z . /EB 

ZF (Z . EQ .2) GO  TO  ZOO 

ZD-SDZ*DPSZ 

ZC-ZD 

EGZI-DCMPLX (DCOS (ZC) , DSZN (ZC) ) 

RR1-R1+ZD-TD1 

RR2-R2+ZD-TD2 

CALL  CZSZ (COSZ, CZH.SZNZ, RRZ) 

CALL  CZSZ(COS2,CZK,SZN2,RR2) 

EXPB-DCMPLX (COS2-COS1, SZNZ-SIN2) 

RR1-R1-ZD+TD1 
RR2— R2-ZD+TD2 

CALL  CZSZ (COSZ, CZN,SZRZ,RR1) 

CALL  CZSZ  (COS2,CZII>  SIN2,RR2) 

EXPA-DCMPLX  (COS2  -COS  Z ,  SZNZ -SZM2 ) 

F  (1,  Z)  -2 . *SGDS* < .0D0, 1 -DO) ‘EXPA/EGZZ 
F (Z , 2) -2 . *SGDS* ( . 0D0 , Z .00) *EXPB«SGZZ 
ZOO  SZ-S2 

CST— ETA/ (Z6 . *PZ*SGDS*SGDT) 

PZZ-CST*  (  (  F  (Z,  Z)  +E (2,2)  «ES2-E  (1, 2)  /ES2)  *ET2 
A+  (-F  (Z,  2)  -E  (2,  Z)  *ES2+E  (Z,  Z)  /ES2)  /ET2) 

RETURN 

ZZO  ZF (CPSZ .LT. 0 . ) GO  TO  Z20 
TA-TZ 
TB-T2 
GO  TO  Z30 
Z20  TA— TZ 
TB— T2 
SGDT— SGDT 
Z30  SZ-SZ 
CZ— CDS 

PZZ— ( . 0D0, .0D0) 

DO  Z50  Z-Z,2 
TJ-TA 

DO  Z40  J-Z,  2 
DZ-TJ-SZ 

R-DSQRT (DSQ+DZ*DZ) 

X-R+DZ 

ZF  (DZ . LT .  . 0)  X-DSQ/  <R-DZ) 

CALL  CZSZ (COSZ, CZM, SZMZ, X) 

EP-DCMPLX (COSZ, -SZNZ) 

X-R-DZ 

ZF  (DZ .  CT . .  0)  X-DSQ/  (R+DZ) 

CALL  CZSZ (COSZ, CZM, SZMZ, X) 

EM- DCMFLX (COSZ , -SZMZ) 

ZF (J.EQ.2)GO  TO  140 

SP1-EP 
EMZ-EM 
Z40  TJ-TB 
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X-TB-SI 

EGDZ-DCMPLX  (DCOS  (X) ,  DSIN  (X)  ) 

Pll-Pll+CI* < (EP-XP1) *EGDZ+ (EM -EMI) /EGDZ) 
CI-1. 

150  SX-S2 

P11-BTA*P11/ (8 . *PI*SGDS*SGDT) 

RETURN 

END 

C 

C 
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COMPUTER  PROGRAM  RICHMOND4  FOR  THE  FAR-ZONE  FIELD  OF  A 
MONOPOLE  ELEMENT  ON  A  DISK  GROUND  PLANE  ABOVE  FLAT  EARTH 


B-l 


COMPUTER  PROGRAM  RICHMOND* 
MONOPOLE  ANTENNA  ON  CIRCULAR  DISK 
OVER  FLAT  EARTH 

(IMPEDANCE,  GAIN,  AND  FAR-FIELD  PATTERNS) 


by 

Jack  H.  Richmond 
February  16, 1990 
Revised  September  21,  1990 


INTRODUCTION  1 

Appendix  I  presents  the  computer  program  RICHMOND4.  This  FOR¬ 
TRAN  program  calculates  the  current  distribution,  terminal  impedance, 
and  directive  gain  G(6)  of  a  base-fed  monopole  antenna  mounted  at  the 
center  of  a  circular  disk  over  the  flat  lossy  earth.  The  detailed  theory  be¬ 
hind  this  moment-method  solution  is  presented  in  the  following  paper: 

(J.  H.  Richmond,  “Monopole  Antenna  on  Circular  Disk  over  Flat  Earth,” 
IEEE  Transactions,  Vol.  AP-33,  pp.  633-637,  June  1985.) 

To  assist  the  user,  comment  statements  have  been  inserted  in  the  main 
program  and  the  subroutines.  Only  a  few  additional  brief  comments  will 
be  required  in  this  Introduction. 

RICHMOND4  performs  all  calculations  with  double  precision.  In  this 
program  the  notation  corresponds  closely  with  the  notation  in  the  above 
paper,  with  one  exception:  In  the  paper  z0  denotes  the  height  of  the  circular 
disk  above  the  surface  of  the  earth,  whereas  in  the  program  HDL  denotes 
Zo! A.  (The  wavelength  in  free  space  is  denoted  by  A  or  WAVM.) 


’Appreciation  U  expressed  to  The  MITRE  Corporation  for  sponsoring  this  report . 
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RICHM0ND4  requires  all  the  subroutines  used  by  RICHMOND3.  In 
addition,  RICHM0ND4  requires  the  following  additional  subroutines  which 
are  listed  after  the  main  program  in  Appendix  I:  EDISKl,  EDISK2,  GAINl, 
and  GAIN2. 

Subroutines  EDISKl  and  GAINl  calculate  the  far-zone  electric  field 
intensity  E$(0)  for  the  monopole  on  the  circular  disk  in  free  space,  with 
0  <  6  <  w-  In  these  calculations,  the  factor  exp(-jkr)/(kr)  is  suppressed. 
For  the  monopole  on  a  circular  disk  over  the  flat  earth,  EDISK2  and  GAIN2 
calculate  Eg(6 )  with  0  <  8  <  tt/2.  In  GAINl  and  GAIN2,  ET  denotes  the 
quantity  kr  exp(jkr)  F«(r,fl),  which  may  be  called  “the  normalized  far- 
zone  electric  field  intensity”  corresponding  to  the  space  wave. 

GAINl  and  GAIN2  each  makes  two  passes  through  the  range  of  angles 
6.  On  the  first  pass  (M  =  1),  the  time-average  radiated  power  Pr  is  calcu¬ 
lated  via  numerical  integration  using  appropriately  small  increments  DTH 
in  the  angle  0.  On  the  second  pass  (AT  =  2),  the  directive  gain  D(0)  is 
calculated  and  stored  using  the  angular  increments  DTHD  specified  in  the 
input  data. 

Tables  I  and  II  show  numerical  results  (with  RICHM0ND4)  for  cir¬ 
cular  disks  with  radii  ka  —  1.5  and  3.0  respectively,  where  k  denotes  the 
wavenumber  in  free  space.  For  the  monopole  on  a  disk  in  free  space,  the  ra¬ 
diation  resistance  R(RAD)  and  the  directive  gain  GAIN  agree  closely  with 
Tables  A2-6  and  A2-12  in  the  following: 

(M.  M.  Weiner  et  al.,  “Monopole  Elements  on  Circular  Ground  Planes,” 
Artech  House,  1987.) 

In  addition,  Tables  I  and  II  list  the  radiation  resistance  and  gain  for  the 
monopole  on  a  disk  on  the  fiat  lossy  earth,  as  well  as  the  antenna  terminal 
impedance  in  free  space  and  on  the  flat  earth.  Let  us  define  the  radiation 
efficiency  to  be  the  ratio  between  the  power  radiated  (via  the  “space  wave” 
into  the  free-space  region)  and  the  power  input  at  the  antenna  terminals. 
This  radiation  efficiency,  then,  is  equal  to  the  ratio  of  the  radiation  resis¬ 
tance  and  the  real  part  of  the  antenna  terminal  impedance.  In  Tables  I 
and  II  the  radiation  efficiency  is  100%  for  the  antenna  in  free  space.  For 
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the  antenna  on  flat  earth  the  radiation  efficiency  is  25.1%  with  the  smaller 
disk,  and  46.6%  with  the  larger  disk. 

RICHMOND4  has  been  tested  with  several  of  the  examples  in  the  pub¬ 
lished  paper  cited  earlier  (Richmond,  IEEE,  1985),  with  excellent  agree¬ 
ment  on  the  antenna  terminal  impedance.  In  addition,  the  RICHMOND4 
calculations  converge  properly  as  the  number  of  segments  (NSD  and  NSW) 
increases. 

In  the  calculation  of  the  normalized  far-zone  electric  field  intensity  ET, 
GAIN1  and  GAIN2  include  only  the  “space  wave”  component  of  the  field. 
The  effects  of  the  round  earth  and  the  ionosphere  are  not  included.  Even 
with  flat  earth,  as  the  observer  approaches  the  air-earth  interface,  the 
“ground  wave”  or  “surface  wave”  field  may  become  significant,  but  is  not 

included  in  ET. 


TABLE  I.  Numerical  Results  with  ka  =  1.5 


DOUBLE  PRECISION 
MONOPOLE  ON  CIRCULAR  DISK 


NSD  NSW 

A L 

CK  CL 

HL  HDL 

8  2 

l.E-6 

1.5  0.2387 

0.25  0.0 

R(RAD) 

R 

X 

22.0733 

22.0710 

13.7758 

(IN  FREE  SPACE) 

THETA 

GAIN 

(IN  FREE  SPACE) 

0.0000 

0.0000 

10.0000 

0.0582 

20.0000 

0.2202 

30.0000 

0.4521 

40.0000 

0.7103 

50.0000 

0.9540 

60.0000 

1.1543 

70.0000 

1.2967 

80.0000 

1.3781 

90.0000 

1 .4008 

100.0000 

1.3674 

110.0000 

1.2775 

120.0000 

1.1306 

130.0000 

0.9303 

140.0000 

0.6906 

150.0000 

0.4387 

160.0000 

0.2135 

170.0000 

0.0564 

180.0000 

0.0000 

THETA 

GAIN 

(ON  FLAT  EARTH) 

0.0000 

0.0000 

10.0000 

0.1550 

20.0000 

0.6084 

30.0000 

1.3146 

40.0000 

2.1628 

50.0000 

2.9360 

60.0000 

3.2924 

70.0000 

2.8331 

80.0000 

1.3814 

90.0000 

0.0000 

NSD  NSW 

AL 

BAR  CL 

ER 

FMC  HDL  HL 

8  2 

l.E-6 

3.0  0.2387 

4. 

10.  0.0  0.25 

R (RAD) 

R 

X 

7.4736 

29.7531 

12.3685 

(ON  FLAT  EARTH) 

/ 
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SIG 

.001 


TABLE  II.  Numerical  Results  with  ka  =  3 


DOUBLE  PRECISION 
MONOPOLE  ON  CIRCULAR  DISK 


NSD  NSW 

AL 

8  2 

1  .E-6 

R (RAD) 

R 

40.0556 

40.0120 

THETA 

GAIN 

0.000 0 

0.0000 

10.0000 

0.2238 

20.0000 

0.7637 

30.0000 

1.3305 

40.0000 

1.6823 

50.0000 

1.7417 

60.0000 

1.5725 

70.0000 

1.2923 

80.0000 

1.0069 

90.0000 

0.7868 

100.0000 

0.6695 

110.0000 

0.6649 

120.0000 

0.7499 

130.0000 

0.8597 

140.0000 

0.8925 

150.0000 

0.7578 

160.0000 

0.4593 

170.0000 

0.1391 

180.0000 

0.0000 

THETA 

GAIN 

0.0000 

0.0000 

10.0000 

0.2630 

20.0000 

0.9524 

30.0000 

1.8271 

40.0000 

2.6187 

50.0000 

3.1019 

60.0000 

3.0919 

70.0000 

2.4352 

80.0000 

1.1249 

90.0000 

0.0000 

NSD  NSW 

AL 

8  2 

1  .E-6 

R(RAD) 

R 

18.2413 

39.1551 

c:t  CL 
3.  0.4775 

X 

33.0500 


BAR  CL  ER 
3.0  .4775  4. 

X 

27.5684 


HL  HDL 

0.2500  0.0 


(IN  FREE  SPACE) 
(IN  FREE  SPACE) 


(ON  FLAT  EARTH) 


FMC  -  HDL  HL 
10.  0.0  .25 


(ON  FLAT  EARTH) 


SIG 

.001 
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APPENDIX  1.  R.ICHMOND4  and  the  Subroutines 


C  RICHMOND 4 

C  MONOPOLE  AT  CENTER  OF  CIRCULAR  DISK  OVER  FLAT  EARTH.  RICKMCH14.1 

C  DOUBLE  PRECISION. 

C  CURRENT  DISTRIBUTION,  IMPEDANCE,  AND  TAR-FIELD  PATTERN. 

C  LINK:  BES10,C1SI,CMINV,DZ11,DZDD,DZWD,DZWW,EDISK1,EDISK2,EXPJ, 

C  GAIN! ,  GA1N2 , GRILL, QDD, QDM, QMM, SKEW , SKEWS , SKEWT , SFART , ZSDM, ZSMM 

IMPLICIT  R£AL*6  (A-H),  (P-Z) 

COMPLEX*  16  CJ  (30)  ,VJ(30)  ,ZJ(30)  ,VIJ(30,30)  ,ZIJ(30,30) 

COMPLEX* 16  Yll , DET, EC, Dll , D12,D21,D22,DZ1J, DV1 ,  Wll 
COMPLEX* 16  P11,P12,P21,P22, ZDD, ZDM, ZMD, ZMM, Zll , ZD, Z22 , Z12, Z21 
DIMENSION  FB  (500)  ,G(182)  ,LLL(30)  ,MMM(30) 

DATA  E0,U0/8 .85418? 33677D- 12,1 .  256633061 44D-6/ 

DATA  ETA,  PI,  TP/ 376. 730366239, 3. 14159265359, 6. 28318530*718/ 

DATA  ICC, IFB/30, 500/ 

1  FORMAT (IX, 215, BF10. 4) 

2  FORMAT (IX, 8F10 . 4) 

5  FORMAT (1H0) 

AL  -  RADIUS  OF  WIRE  IN  WAVELENGTHS. 

BAR  *  RADIUS  RATIO  FOR  COAXIAL  FEED  CABLE. 

CL  -  RADIUS  OF  CIRCULAR  DISK  IN  WAVELENGTHS  =  EPSLN/TP. 

DTHD  -  INCREMENT  IN  FAR-FIELD  ANGLE  THETA  (DEGREES) . 

ER  -  RELATIVE  PERMITTIVITY  OF  EARTH. 

FMC  «  FREQUENCY  IN  MEGAHERTZ. 

HL  -  LENGTH  OF  MONOPOLE  IN  WAVELENGTHS. 

HDL  -  HEIGHT  OF  DISK  ABOVE  THE  EARTH  IN  WAVELENGTHS. 

NSD  »  NUMBER  OF  SEGMENTS  ON  THE  DISK. 

NSW  «  NUMBER  OF  SEGMENTS  ON  THE  WIRE. 

SIG  «  CONDUCTIVITY  OF  EARTH,  MHO/M. 

SET  DTHD  -  NEGATIVE  TO  SKIP  THE  GAIN  CALCULATIONS. 

SET  HDL  -  NEGATIVE  FOR  MONOPOLE-DISK  IN  FREE  SPACE, 

OR  HDL  -  POS.  FOR  FREE  SPACE  +  FLAT  EARTH. 

SET  IWCJ  -  1  TO  WRITE  OUT  THE  CURRENTS  CJ(N), 

OR  IWCJ  «  0  TO  SUPPRESS  WRITEOUT . 

TL  *  l.D-5  FOR  EPSLN  GREATER  THAN  OR  EQUAL  0.25, 

«  AL.D-4  FOR  EPSLN  LESS  THAN  0.25. 

AL-l.D-6 
BAR=3  . 

CL* 3  .  /TP 
DTHD* 10 . 

ER=4  . 

FMC-10 . 

HL* . 25 
HDL-l.D-5 

iwcj-0 

NSD«8 
NSW-2 
SIG* . 001 
TL-1 .D-5 
WAVM-300  .  /FMC 
IWZ-0 
NPH-6 

NEQ-NSD+NSW-1 
IF (NEQ. GT. ICC) GO  TO  400 
AK-TP*AL 
CK-TP*CL 
HK*TP‘HL 
HDK-TP*HDL 
TK-TP*TL 
OMEG-TP  *  FMC  *  1 .  D  6 
EC-DCMPLX (ER, -SIG/ (OKEG*EO> ) 

DKD-(CK-AK) /NSD 
DKW-HK/NSW 
RH2-AK+DKD 

IF  (RH2.LT.  BAR*  AK)  GO  TO  400 
TDKD-2 . *DKD 
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CDKD-DCOS (DKD) 

SDKD-DSIN (DKD) 

CDK-DCOS (DKW) 

SDK-DSIN (DKW) 

MXX-NSW-1 

NA-NSD+1 

CALL  QMM  (AX,  DKD,  DKW,  CDKD,  SDKD,  CDK,  SDK,  TX,  IWZ, NPH,  211) 

ZIJ(1,1)-Z11 

IF (MSD . LE . 1 ) GO  TO  100 

51- AK 

DO  60  J-2,NSD 

52- S1+DKD 

53- S1+TDKD 
Tl-AK 

DO  50  1-2, J 
T2-T1+DKD 
‘  T3— Tl+TDKD 

CALL  QDD  (CDKD,  SDKD,  SI,  S3,T1,T3,  TK,  IWZ, NPH,  222) 

ZIJ(I, J)-Z22 
50  T1-T1+DKD 

CALL  QDM  (AK, DKD,  DKW,  CDKD,  SDKD,  SDK,  SI,  S3,  TK,  IWZ, NPH,  Z12) 
ZIJ (1, J) — Z12 
60  SI— Sl+DKD 
100  IF (NSW . LE . 1 ) GO  TO  200 

CALL  SPART  (AK, DKD, DKW, MAX,  IWZ,  ZJ, CJ) 

L— 0 

DO  160  I— NA, NEQ 
DO  150  J— I,NEQ 
K-J-I+l 

150  ZIJ ( I , J) — Z J (K) 

L-L+l 

ZIJ (1,1) —CJ (L) 

160  CONTINUE 

178  IF (NSD . LE . 1 ) GO  TO  200 
Z2-.0 

DO  190  J— NA, NEQ 
Z2— Z2+DKW 
SI— Z2-DKW 
S3— Z2+DKW 
RH2-AK 

DO  180  1-2, NSD 
RH2  — RH2  +DKD 
T1-RH2-DKD 
T3-RH2+DKD 

CALL  SKEWT  (AK,  SI,  S3,  Tl, T3,  CDK,  SDK,  CDKD,  SDKD,  IWZ,  Z12) 
180  ZIJ (I, J) — Z12 
190  CONTINUE 

200  CALL  GRILL  (AK,  BAR,  DKD,  DKW, NEQ,  NSD,  NSW,  TK,VJ) 

DO  210  I— 1 , NEQ 
DO  210  J-I,NEQ 

C  WRITE  (17,1)  I,  J,  ZIJ  (I,  J) 

ZIJ  (J,  I) -ZIJ  (I,  J) 

210  VIJ (I, J) -ZIJ (I, J) 

C  WRITE  (6,1)  NSD,  NSW,  AL,CK,  CL,  HL,HDL 

C  WRITE  (17, 1)NSD,NSW,AL,CK,CL,HL,HDL 

C  WRITE  (6,  5) 

C  WRITE  (17,5) 

CALL  CMINV(CJ,VJ,  ZIJ,  ICC,  IWCJ,  1,LLL,MMM,KEQ,DET) 

Til— CJ  (1) 

zii-i./ni 

c 

C  CALCULATE  DIRECTIVE  GAIN  G(N)  IN  FREE  SPACE. 

C 

RR1-.0 

C  IF (DTHD.LE. .0) GO  TO  212 

c  .  GAINA  *  DIRECTIVITY  IN  FREE  SPACE. 
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C  THA  -  ANGLE  OF  MAXIMUM  GAIN  IN  rREE  SPACE. 

CALL  GAIN1  (AK,  CK,  CJ,  DTHD,  G,  GAINA,  HK,  NSD, 

2  NSW, NTH, PR, THA, WAVM) 

AJ-CDABS (CJ(1) ) 

RR1-PR/  <AJ*AJ) 

212  CONTINUE 

WRITE  ( 6 , 2 )  HDL ,  GAINA,  RR1 ,  Z 1 1 
WRITE  (15, 2)  HDL,  GAINA, RR1,Z11 
WRITE  (6,  5) 

WRITE  (17,  5) 

IF (DTHD . LE . . 0 ) GO  TO  222 
DO  220  N-1,NTH 
TH-DTHD* (N-l) 

WRITE  (6,2)  TH,  G  (N) 

WRITE  (17,2)  TH, G  (N) 

220  CONTINUE 
WRITE (6,5) 

WRITE  (17, 5) 

222  CONTINUE 

CALCULATE  DIRECTIVE  GAIN  G (N)  OVER  FLAT  EARTH. 

IF (HDL . LT . . 0 ) GO  TO  350 
DO  400  NHD-1,4 
HDL-NHD 
HDK-TP*HDL 

C  DELETE  STATEMENT  230  UNLESS  THE  CURRENT  DISTRIBUTION  IS  TO  BE 
C  APPROXIMATED  BY  THE  CUR.  DIST.  FOR  ANTENNA  IN  FREE  SPACE. 

C  230  IF (NHD . GT . 0 ) GO  TO  316 

CALL  DZ11  (AX,BAR, DKD,DKN,EC,FB,  HDK,TK,  IFB,D11,DV1) 

ZIJ (1, 1) “Dll 
IF (NSD . LE . 1 ) GO  TO  265 
S2-AK+DKD 
DO  260  J«2,NSD 
T2-AK+DKD 
112-1 

DO  250  1-2, J 

CALL  DZDD  (AK,DBET,DKD,DKW,EC,FB,HDK,  S2,T2,TK 
2, IFB, 112, KMX, D12, D22) 

IF (I .EQ . 2) P12-D12 
ZIJ (I, J)-D22 
112-2 

250  T2-T2+DKD 

ZIJ(1, J)-P12 
260  S2— S2+DKD 
265  IF(NSW.LE.l)GO  TO  276 
DO  276  K— 1,MAX 

CALL  DZWW  (AK,  DKD,DKW,EC,HDK,K,TK,  ZJ,DZ1J) 

J-NA+K-1 
ZIJ(1, J)-DZ1J 
L-l 

DO  270  I— NA, J 
ZIJ(I,J)  -ZJ  (L) 

270  L-L+l 
276  CONTINUE 

276  IF (NSD . LE . 1 ) GO  TO  300 
Z2-.0 

DO  290  J-NA, NEQ 
Z2-Z2+DKW 

CALL  DZND  (AK,DKD,DKW,EC,HDK,NSD,TK,  Z2,  ZJ) 

DO  280  1-2, NSD 
280  ZIJ(I,J)-ZJ(I) 

290  CONTINUE 
300  DO  310  1-1 , NEQ 
DO  308  J-I , NEQ 
Z12-VTJ(I, J) 
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D12-ZIJ  (I,  J) 

C  WRITE (17,1) I, J,Z12,D12 

ZIJ (I, J) -Z12+D12 
308  CONTINUE 
310  CONTINUE 
C  WRITE  (17, 5) 

CALL  GRILL  (AK,BAR,DKD,DKW,NEQ, NSD, NSW, TK,VJ) 
VJ(1)-VJ(1)+DV1 
DO  315  I-1,NEQ 
DO  312  J-I,NEQ 
312  ZIJ  (J,  I)  “ZIJ  (I,  J) 

315  CONTINUE 

CALL  CMINV  <CJ,  VJ,  ZIJ,  ICC,  IWCJ,  1,  LLL, MMM,NEQ,  DET) 

316  RR2-.0 

C  IF (DTHD.LE. . 0) GO  TO  322 

C  GAINS  -  DIRECTIVITY  OVER  FLAT  EARTH. 

C  THB  -  ANGLE  OF  MAXIMUM  GAIN  OVER  FLAT  EARTH. 

CALL  GAIN2  (AK,CK,CJ,DTHD,EC,G,GAINB,HDK,HK, 

2  NSD , NSW , NTH, PR,  THB , WAVM) 

AJ-CDABS (CJ(1) ) 

RR2-PR/  (AJ*AJ) 

IF (DTHD . LE . 0 . ) GO  TO  322 
DO  320  N-1,NTH 
TH-DTHD* (N-l) 

WRITE  (6, 2)  TH,G  (N) 

WRITE  (17,2)  TH,  G  (N) 

320  CONTINUE 
WRITE (6, 5) 

WRITE (17,5) 

322  Yll-CJ(l) 

Wll-l./Xll 

DBB-10 . *ALOG10 (GAINB) 

WRITE ( 6 , 2 ) HDL , DBB , THB , RR2 , W1 1 
WRITE (1 6 , 2 ) HDL, DBB , THB , RR2 , W1 1 
C  WRITE  (17,1)  NSD ,  NSW,  AL,  BAR,  CL,  ER,  FMC,  HDL,  HL,  SIG 

C  WRITE (17,5) 

350  CONTINUE 
400  CONTINUE 

WRITE (6, 5) 

WRITE  (16,  5) 

DBA-1 0.*ALOG10 (GAINA) 

WRITE (6, 2) DBA, THA, RR1,  Zll 
WRITE  { 1 6 , 2 )  DBA,  THA,  RR1 ,  Zll 
500  CALL  EXIT 
END 
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SUBROUTINE  GAIN1  (AX,  CK,CJ,  DTHD, 6, GAINA, BK, NSD, NSW, 

2  NTH, PR, THA, WAVM) 

SEPTEMBER  21,  1990. 

CALCULATE  DIRECTIVE  GAIN  G (N)  FOR  MONOPOLE  ON  DISK  IN  FREE  SPACE. 
ALSO  PR  «  TIME-AVERAGE  POWER  RADIATED. 

GAINA  -  DIRECTIVITY  IN  FREE  SPACE. 

THA  -  ANGLE  OF  MAXIMUM  GAIN  IN  FREE  SPACE. 

IMPLICIT  REAL* 8  (A-H) ,  (P-Z) 

COMPLEX*  16  CJ(1)  ,CJ1,CJ2,CQ,CQS,EK1,EK2,EK3,ET,ETH,ETHD 
DIMENSION  G (1) 

DATA  ETA, PI, TP/376 .730366239, 3 . 14159265359, 6.28318530718/ 
BET-TP/WAVM 
DKD- (CK-AK) /NSD 
SDKD-DSIN (DKD) 

CDKD-DCOS (DKD) 

DKW-HK/NSW 
CDK-DCOS (DRW) 

SDK-DSIN(DKW) 

NEQ-NSD+NSW-1 
DK-HK 

IF (CK . GT . HK) DK-CK 
NS«10.*DK/PI 
NS-2* (NS/2) 

IF (NS.LT.20)NS-20 
NTH-NS+1 
DTH-PI/NS 

CQ-DCMPLX  ( .  0D0 , BET*ETA/  (4 .  *PI*SDK)  ) 

ZO-HDK 
ZO-.O 

LS-10.*DKD/PI 
IF(LS.LT.4)LS-4 
LS-2* (LS/2) 

PR-.O 

GAINA- . 0 
DO  250  M-1,2 

IF  M-l,  USE  SIMPSON'S  RULE  INTEGRATION  TO 
CALCULATE  PR  -  POWER  RADIATED. 

IF  H-2,  CALCULATE  THE  DIRECTIVE  GAIN  G  (N)  . 

SGN— 1 . 

DO  200  NT-1, NTH 
IF (NT . EQ . 1 ) GO  TO  200 
IF (NT. EQ. NTH) GO  TO  200 
ET- ( . 0D0, .0D0) 

WF-3.+SGN 

C  SELECT  THE  FAR-FIELD  ANGLE  TH  -  THETA. 

TH-DTH* (NT-1) 

CTH-DCOS (TH) 

STH-DSIN (TH) 

CQS-CQ/STH 

C  CALCULATE  FAR- ZONE  FIELD  OF  MODE  fl  IN  FREE  SPACE. 

Z1-Z0 
Z2-Z1+D1CW 
ARC— Z2*CTH 

EK2— DCMPLX  (DCOS  (ARG)  ,  DSIN  (ARG) ) 

ARC— Z1*CTH 

EK1-DCMPLX  (CDK,  CTH*SDK)  *  DCMPLX  (DCOS  (ARG)  ,  DSIN  (ARG)  ) 

ETH-CQS  *CJ (1 ) * (EK2-EK1 ) 

CJ1— CJ(1) 

CJ2- ( . 0D0 , .0D0) 

RK1-AX 

RK2-RK1+DKD 

CALL  EDISK1  (CJ1,CJ2,CTH,DKD,ETHD,LS,RX1,RX2,SDKD,STH,BET,Z0) 
ET-ET +ETH+ETHD 
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IF (NSD . LE . 1 ) GO  TO  100 

CALCULATE  FAR  FIELD  FROM  DISK  CURRENT  IN  FREE  SPACE. 

DO  60  J— 1,NSD 
RK1«AK+(J-1)*DKD 
RK2-RK1+DKD 
CJl-CJ(J) 

IF ( J . EQ . 1 ) C  Jl- ( . 0D0 , .0D0) 

CJ2-CJ(J+1) 

IF (J.EQ.NSD) CJ2- ( . 0D0, .0D0) 

C  EDISK1  CALCULATES  FIELD  FROM  ANNULAR  ZONE  J  OF  THE  DISK. 

C  RX1  AND  RK2  ARE  INNER  AND  OUTER  RADII  OF  ZONE  J. 

C  CJ1  AND  CJ2  DENOTE  RADIAL  CURRENTS  AT  INNER  AND  OUTER  RADII. 

CALL  EDISK1  (CJ1 ,  CJ2,  CTH,  DKD, ETHD,  LS,RK1,RK2,  SDKD,  STH,BET,  Z0) 
ET-ET+ETHD 
60  CONTINUE 

100  IF (NSW . LE . 1 ) GO  TO  162 

C  CALCULATE  FAR  FIELD  FROM  THE  WIRE  DIPOLE  MODES  IN  FREE  SPACE. 
JB-NSW-1 
DO  160  J-l, JB 
Z2-Z0+J*DKW 
Z1-Z2-DKW 
Z3-Z2+DKW 
L-NSD+J 
ARG-Z1 *CTH 

EK1-DCMPLX  (DCOS  (ARG)  ,  DSIN  (ARC)  ) 

ARG-Z2*CTH 

EK2-DCMPLX  (DCOS  (ARG)  ,DSIN  (ARG)  ) 

ARG-Z3*CTH 

EK3-DCMPLX  (DCOS  (ARG)  ,  DSIN  (ARG)  ) 

ETH-CJ  (L)  *CQS*  (EK1+EK3-2 .  *CDK*EK2) 

ET-ET+ETH 
160  CONTINUE 
162  AT-CDABS(ET) 

ATS-AT**2 

IF (M . EQ . 1 ) PR-PR+WF*ATS*STH 
IF  (M .  EQ .  2  )  G  (NT  )  -ATS 

ZF (ATS . LT . GAINA) GO  TO  200 
GAINA— ATS 
THA«57.29578*TH 
200  SGN— SGN 

IF  (M.EQ . 1) PR-TP*PR*DTH/  (3 .  *ETA*BET*BET) 

DTH- . 0 1 7 4 5 32 9 *DTHD 

IF (DTHD . GT . 0 . ) GO  TO  248 
NS— 1 
NTH— 1 
GO  TO  250 

248  NS— 180 . /DTHD 

NTH-NS+1 
250  CONTINUE 
G(l)-.0 
G  (NTH)  — .  0 

c  NORMALIZE  THE  DIRECTIVE  GAIN  G(N)  . 

CST-4 . *PI/ (ETA*BET*BET*PR) 

GAINA— CST*GAINA 
DO  300  N-l , NTH 
GN— CST*G (N) 

C(N)-GN 
300  CONTINUE 
RETURN 
END 
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SUBROUTINE  GAIN2  (AX,  CK,  CJ, DTHD,  EC,  G,  GAINB,  HDK,  HK,  GAB2.1 

2  NSD,  NSW, NTH, PR, TBB, WAVM) 

SEPT .  21  1990. 

GAINB  -  DIRECTIVITY  OVER  FIAT  EARTH. 

THB  -  ANGLE  OF  MAXIMUM  GAIN. 

CALCULATE  DIRECTIVE  GAIN  G(N)  FOR  MONOPOLE  ON  DISK  OVER  FLAT  EARTH. 
ALSO  PR  -  TIME-AVERAGE  POWER  RADIATED  INTO  UPPER  HALF-SPACE. 

DOUBLE  PRECISION 

IMPLICIT  REAL*8  (A-H) ,  (P-Z) 

COMPLEX* I 6  CJ<1) ,CJ1,CJ2,CQ,CQS 

COMPLEX* 1 6  EC, EK1 , EK2 , EK3 , ET, ETH, ETHD, ET1 ,  QST, RC 
DIMENSION  G (1) 

DATA  ETA,  PI,  TP/376 . 730366239, 3 . 14159265359, 6.28318530718/ 

.  BET-TP/WAVM 
DKD- (CK-AK) /NSD 
SDKD-DSIN(DKD) 

CDKD-DCOS (DKD) 

DKW-HK/NSW 
CDK-DCOS (DKW) 

SDK-0SIN (DKW) 

NEQ— NSD+NSW- 1 

DK-HK+HDK 

IF (CK . GT . DK) DK-CK 

NS-10.*DK 

NS-2* (NS/2) 

IF (NS.LT.90)NS-90 

NTH-NS+1 

DTH-PI/(2.*NS) 

CQ-DCMPLX ( . 0D0 , BET*ETA/ (4 . *PI*SDK) ) 

ZO-HDK 

LS-10.*DKD/PI 
IF(LS.LT.4)LS-4 
LS-2* (LS/2) 

PR-.O 

GAINS-. 0 
DO  250  M-1,2 

IF  M-l,  USE  SIMPSON'S  RULE  INTEGRATION  TO 
CALCULATE  PR  -  POWER  RADIATED. 

IF  M-2,  CALCULATE  DIRECTIVE  GAIN  G(N). 

SGN--1 . 

DO  200  NT-1, NTH 
IF (NT . EQ . 1 ) GO  TO  200 
ET— ( . 0D0, .0D0) 

WF-3.+SGN 

IF (NT . EQ . NTH) WF-1 . 

SELECT  THE  FAR-HELD  ANGLE  TB  -  THETA. 

TH-DTH* (NT-1) 

CTH-DCOS (TH) 

STH-DSIN (TH) 

CQS-CQ/STH 

QST-CDSQRT (EC-STH*STH) 

RC  -  REFLECTION  COEFFICIENT  AT  AIR-EARTH  INTERFACE. 

RC- (EC*CTH-QST) / (EC*CTH+QST) 

CALCULATE  FAR- ZONE  FIELD  OF  MODE  #1. 

Z1-Z0 

Z2-Z1+DKW 

ARG1“Z1*CTH 

EK1-DCMPLX  (DCOS  (ARG1)  ,  DSIN  (ARG1)  ) 

ET1-DCMPLX (CDK, CTH*SDK) *EK1 
ARG2-Z2 *CTH 

EK2-DCMPLX  (DCOS  (ARG2)  ,D3IN  (ARG2)  } 
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no  onoo  on 


GA2H2.2 


ETH-EX2-ET1 
EKl-DCONJG(EKl) 

ET1-DCMPLX (COX, -CTH*SDK) *EK1 
EK2-DC0NJG(EK2) 

£TH-CQS*CJ (1) * (ETH+RC* (EK2-ET1) ) 

CJ1— CJ<1) 

CJ2- ( . 0D0,  .0D0) 

RK1-AK 
RK2-RK1+DKD 

CALL  EOISK2  (CJ1 ,  C  J2 ,  CTH ,  DKD ,  ETHD ,  LS ,  RC ,  RK1 ,  RX2 ,  SDKD ,  STH , BET ,  Z 0 ) 

ET-ET+ETH+ETHD 

IF (NSD . LE . 1 ) GO  TO  100 

CALCULATE  FAR  FIELD  FROM  DISK  CURRENT  IN  FREE  SPACE. 

DO  60  J-1,NSD 
RK1-AK+ (J-1)*DKD 
RK2-RK1+DKD 
CJ1-CJ ( J) 

IF(J.EQ.1)CJ1-(.ODO, .0D0) 

CJ2«CJ(J+1) 

IF ( J.EQ.NSD) CJ2— ( . 0D0 , .ODO) 

EDISK2  CALCULATES  FIELD  FROM  ANNULAR  ZONE  J  OF  THE  DISK. 

RK1  AND  RK2  ARE  INNER  AND  OUTER  RADII  OF  ZONE  J. 

CJ1  AND  CJ2  DENOTE  RADIAL  CURRENTS  AT  INNER  AND  OUTER  RADII. 

CALL  EDISK2  (CJ1 ,  CJ2 ,  CTH, DKD, ETHD, LS, RC,  RKI , RK2 ,  SDKD,  STH, BET,  Z0) 
ET-ET+ETHD 
60  CONTINUE 

100  IF (NSW . LE . 1 ) GO  TO  162 

CALCULATE  FAR  FIELD  FROM  THE  WIRE  DIPOLE  MODES  IN  FREE  SPACE. 

JB-NSW-1 

DO  160  J-l, JB 

Z2«Z0+J*DKW 

Z1-Z2-DKN 

Z3-Z2+DKW 

L-NSD+J 

ARG-Z1*CTH 

EK1-DCMPLX (DCOS (ARG) , DSIN (ARG)  } 

ARG-Z2*CTH 

EK2-DCMPLX (DCOS (ARG) , DSIN (ARG)  ) 

ARG«Z3*CTH 

EK3-DCMPLX (DCOS (ARG) , DSIN (ARG) ) 

ETH-EK1+EK3-2 . *CDK*EX2 
EKl-DCONJG(EKl) 

EK2-DCONJG (EK2) 

EK3-DCONJG (EK3) 

ETH-CJ  (L)  *CQS*  (ETH+RC*  (EK1+EK3-2  .  *CDK*EK2)  ) 

ET-ET+ETH 
160  CONTINUE 
162  A1>CDABS(ET) 

ATS-AT**2 

IF (M.EQ. 1) PR-PR+NF*ATS*STH 
IF  (M.EQ.2)G (NT) -ATS 

IF (ATS.LT.GAINB) GO  TO  200 
GAINB-ATS 
THB— 57 . 2  957  8  *TH 
200  SGN—  SGN 

IF  (M .  EQ .  1 )  PR-TP *PR*DTH/  ( 3 .  *ETA*BET*BET) 

IF (DTHD . GT . 0 . ) GO  TO  248 

NS-1 

NTH— 1 

GO  TO  250 

248  DTH— .01745329*DTHD 
MS— 90. /DTHD 
NTH-NS+1 
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250  CONTINUE 

G(1)«.0  GA  1X2.3 

C  NORMALIZE  THE  DIRECTIVE  CAIN  C(N)  . 

CST-4 .  *PI/ <ETA*BET*BET*PR) 

GAINB-CST *  GAINB 
DO  300  N-1,NTH 
GN-CST*G (N) 

G(N)«GN 
300  CONTINUE 
RETURN 
END 
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APPENDIX  C 


COMPUTER  PROGRAM  RICHMD6  FOR  THE  INPUT  IMPEDANCE,  CURRENT 
DISTRIBUTION,  AND  FAR-ZONE  FIELD  OF  A  MONOPOLE  ELEMENT  ON  A 

PERFECT  GROUND  PLANE 


C-l 


C  PROGRAM  "RICHMD6 . FOR" 

C***** ************************************************************************ 

c  * 

c  * 

C  THIS  COMPUTER  PROGRAM,  IN  FORTRAN  LANGUAGE,  VAS  WRITTEN  BY  DR.  * 

C  JACK  RICHMOND  OF  OHIO  STATE  UNIVERSITY.  IT  USES  A  SINUSOIDAL-GALBRKIN  * 

C  METHOD  OF  MOMENTS  TO  COMPUTE  THE  INPUT  IMPEDANCE,  CURRENT  DISTRIBUTIONS,  * 
C  AND  ANTENNA  PATTERN  OP  A  MONOPOLE  ELEMENT  OF  LENGTH  h  AND  RADIUS  b  ON  A  * 
C  PERFECT  GROUND  PLANE  OP  INFINITE  EXTENT.  A  DETAILED  DERIVATION  IS  * 

C  PUBLISHED  IN  REFERENCE  1.  * 

C  * 

C  * 

C  REFERENCE:  * 

C  * 

C  1)  J.H.  RICHMOND,  "COMPUTER  PROGRAM  VAIT-SURTEES",  REPORT  * 

C  PREPARED  POR  THE  MITRE  CORPORATION,  29  DEC.  1989.  * 

C  .  * 

C***************************************************************************** 

c  * 

C  THIS  COMPUTER  PROGRAM  REQUIRES  FIVE  INPUTS  WHICH  ARE  ENTERED  PROM  * 

C  AN  INPUT  FILE  NAMED  "RICH6_IN.DAT".  * 

C  * 

C  * 

C  b/A  -  AL  -  MONOPOLE  ELEMENT  RADIUS  IN  WAVELENGTHS  * 

C  h/A  -  HL  -  MONOPOLE  ELEMENT  LENGTH  IN  WAVELENGTHS  * 

C  b./b  -BAR  -  RATIO  OF  OUTER  TO  INNER  CONDUCTOR  RADII  OF  * 

C  THE  COAXIAL  LINE  FEED.  * 

C  <7W  -  CMM  .  CONDUCTIVITY  OF  MONOPOLE  ELEMENT  * 

C  (MEGAMHOS/METER)  * 

C  -  -1  FOR  PERFECTLY  CONDUCTING  MONOPOLE  ELEMENT  * 

C  I FLAG  .  FLAG  FOR  MONOPOLE  ELEMENT  CURRENT  DISTRIBUTIONS  * 

C  -  0,  MONOPOLE  ELEMENT  CURRENTS  COMPUTED  BY  METHOD  OF  * 

C  MOMENTS  * 

C  -  -1,  MONOPOLE  ELEMENT  CURRENTS  WITH  IDEALIZED  * 

C  SINUSOIDAL  DISTRIBUTION  AND  FOR  INFINITE  * 

C  CONDUCTIVITY  OF  THE  MONOPOLE  ELEMENT  * 

C  f  -  FMC  -  FREQUENCY  IN  MEGAHERTZ  * 

C  NOTE:  IF  CMM  .  -1,  IT  IS  NOT  NECESSARY  TO  * 

C  SPECIFY  AN  INPUT  VALUE  FOR  FREQUENCY  * 

C  * 

C  * 

C***************************************************************************** 

c 

c  LINK  CISI ; FRILLS ; TPLZ ; TSPAR j ZSURF 


REAL  G(100) , ANG(IOO) .DIRECDB(IOO) .DIRECRP(IOO) .DIRMAX 

INTEGER  HALFNTB.NTH 

COMPLEX  CQQ,CJI,CJA,CJB,DELZ 

COMPLEX  EC,EGZ,EJH,EJS,ETA2,ETD,ETH 

COMPLEX  GM,GP,Q1,Q2,Q3,RC 

COMPLEX  Y11,Z11.ZFIN,ZINF,ZS,ZSG 

COMPLEX  CJ(99),GI(99),U(99),VJ(99),W(99),ZJ(99) 

OPEN (UNIT- 4 , FILE- ' RICH 6  IN . DAT ' , STATUS- 'OLD' , READONLY) 

OPEN (UNIT-10, FILE- 'RICB5  OUT.DAT' , STATUS- 'NEW') 

DATA  ETA, PI , TP/376 . 730365239 , 3 . 14139265359 , 6 . 28318530718/ 

DATA  BO , UO/8 . 85418533677B-12 , 1 . 25663706144E-6/ 

DATA  P2 , B/l . 57079632679 ,2.718281828/ 

DATA  I DM/99/ 

C*********************** ****************************************** *********** 

C  WAVM  -  WAVELENGTH  IN  FREE -SPACE  (METERS).  * 

C  NOTE:  NOT  USED  IF  CMM  -  -1.  * 
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c**************************************************************************** 

2  F0RMAT(7X,8F11.5) 

3  F0RMAT(8X,8F16. 10) 

4  F0RMAT(4X,8F11 .5) 

7  FORMAT(SX) 

521  FORKAT(IX.A) 

C**************************************************************************** 

C  INPUTS  * 

C**************************************************************************** 
R£AD( 4 , * ) AL , HL , BAR , CNM , IFLAG 
IF(CMM.EQ.-1)G0T0  8 
READ(4,*)  PNC 
WAVM-300. /FMC 

8  CONTINUE 
IVR.1 

C*********************************************************** ***************** 

C  INPUT  WRITE  STATEMENTS  * 

C**************************************************************************** 

CALL  HEADING2 ( AL , HL , BAR , CMM , IFLAG ) 

AK-TP*AL 

HK-TP*HL 

SK-2.*HK 

EJS-CMPLX(COS(SK) ,SIN(SK) ) 

CHK-COS(HK) 

SHK-SIN(HK) 

EJH-CMPLX(CHK, SHK) 

NS-15.*BL 

IF(NS.LT.6)NS-6 

IF(NS.GT.IDM)NS-IDM 

IG-NS/2 

IF( IFLAG. EQ.O)GOTO  43 

IG-1 

CNN— 1. 

43  CONTINUE 
NS-2*IG 
N-NS-1 
ZS-(.0,.0) 

IF(CMM.GT.O. )CALL  ZSURF(AK,CMM,FNC,ZS) 

C*****t*********************************************^***tH********ti******** 

C  ZS  -  SURFACE  IMPEDANCE  OF  MONOPOLE  WIRE.  * 

C**  *********************************************  ****************************** 

DK-HK/IG 

CDK-COS(DK) 

SDK-SIN(DK) 

C***************************************************************************** 

ZJ(N)  .  FIRST  ROW  OF  IMPEDANCE  MATRIX  FOR  DIPOLE  IN  FREE  SPACE,  * 

USING  GALERKIN'S  METHOD  VITB  OVERLAPPING  SINUSOIDAL  BASIS  FUNCTIONS.  * 
(DIPOLE  LENGTH  -  TWICE  THE  MONOPOLE  LENGTH,  USING  IMAGE  THEORY.)  * 

C************************** ******************************* ******************** 
CALL  TSPAR(AK,DK,N,ZJ) 

IF(CMM.LT.O. )GO  TO  52 
GK-2.*(DR-CDR*SDR) 

GUSDK-DK*CDK 
FH-4 . *PI*AX*SDK*SDR 
ZJ(1)-ZJ(I)*ZS*QC/FH 
ZJ(2)-ZJ(2)+ZS*GL/FH 
52  112-1 


C  SET  UP  THE  MOMENT-METHOD  VOLTAGE  COLUMN  VJ(M)  FOR  DIPOLE  IN  FREE  SPACE.* 

C************************************************* **************************** 


CALL  FRILLS(AK,BAR,DK,N,CJ) 

DO  60  I-1,N 
K«1+IABS(IG-I) 

60  VJ(I)-CJ(K) 

C ***************************************************************************** 

C  SOLVE  THE  SIMULTANEOUS  LINEAR  EQUATIONS  TO  DETERMINE  THE  CURRENTS  CJ(N>* 

C  ON  THE  WIRE  DIPOLE  IN  FREE  SPACE.  (THE  IMPEDANCE  MATRIX  IS  TOEPLITZ.)  * 

q*****  *******************************************************************  ****** 
CALL  PRCURR 

CALL  TPLZ(CJ,U,VJ,W,ZJ,IER,IWR,I12,N) 

VRITE(10f 521)CHAR(12) 

C***************************************************************************** 

C  ZINF  -  IMPEDANCE  OF  MONOPOLE  ANTENNA  OVER  INFINITE  GROUND  PLANE.  * 

C***************************************************************************** 

CALL  HEADING2(AL,BL,BAR,CMM,IFLAG) 

ZINP-.5/CJ(IG) 

RIN-REAL(ZINF) 

XIN-AIMAG(ZINF) 

CALL  PRRES(RRAD,RIN,XIN) 

Yll-CJ(IG) 

PD-.O 

PIN-REAL (Yll)- 
PRAD-PIN-PD 
EFF-100.*PRAD/PIN 
DTH-2. 

NTH-1.+90./DTH 

DO  I-l.NTH 

ANG(I)»(I-1)*DTH 

TH-ANG(I) 

CALL  GAIN(CDK,CJ,DGAIN,DK,U,HK,N,PRAD,SDK,TH,VJ,WAVH) 

G(I)-DGAIN*2 
END  DO 

DO  I-(NTH+1),91 
G(I)-0 

ANG(I)-(I-1)*DTH 
END  DO 

C  DO  1-1 , NTH 

C  WRITE(10, 2)ANG(I) ,G(I) , ANG(I+NTH) ,G(I+NTH) 

C  END  DO 

C****************************************************************************** 

C  OUTPUT  STATEMENTS  FOR  FAR-FIELD  PATTERN  HEADINGS  (WITH  EARTH)  * 

Q****************************************************************************** 

VRITE(10,850) 

850  F0RMAT(//1X, 

A  ' ELEVATION ' , 5X , ' DIRECTIVE ' , 6X , ' DIRECTIVE ' , 5X , ' RELATIVE  ' ,  12X , 

A  'ELEVATION',  5X,  'DIRECTIVE'  ,6X,  'DIRECTIVE',  5X,  'RELATIVE'  ) 

VRITEC10, 852) 

852  FORMAT (3X, 'ANGLE' ,9X, 'GAIN' , 11X, 'GAIN' ,9X, 'POWER' , 

A  16X, 'ANGLE' ,9X, 'GAIN' ,11X, 'GAIN' ,9X, 'POWER'  ) 

WRITE(10,854) 

854  FORHAT(3X,' (DEG)', 7X,' (NUMERIC) ' ,8X, ' (DBI) ' ,9X, '(DB)' , 

A  16X, '(DEG)', 7X, '(NUMERIC)', 8X, '(DBI)', 9X,'(DB)',/5X  ) 

C I***************************************************************************** 

c****************************************************************************** 
DIRMAX-0.0 
DO  N-l.NTB 
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IF(G(N).NE.0)DIRECDB(N)«10.0*L0G10(G(N)) 

IF ( DIRMAX . LT . DIRECDB (  N ) ) DIRMAX- DIRECDB ( N ) 

END  DO 
DO  N-l.NTH 

DIRECRP ( N ) -DIRECDB ( N ) -DIRMAX 
END  DO 

C  OUTPUT  STATEMENTS  FOR  FAR  FIELD-PATTERN  (VITH  EARTH)  * 

C****************************************************************************** 


HALFNTH  -  NTH/ 2 

DO  N-l, HALFNTH 

IF(  G(N).NE.O)  THEN 

VRITE(10,750)ANG(N),G(N),DIRECDB(N),DIRECRP(N), 

&  ANG  (  N+ HALFNTH  ) ,  G  (  N-f HALFNTH ) ,  DIRECDB  (N-f HALFNTH  ) , 

i,  DIRECRP(N+HALFNTH) 

750  FORMAT(3X,F4.0,2F16.5,F13.5,14X,F4.0,2F16.5,F13.5) 

ELSE 

VRITE(10,604)ANG(N) ,ANG(N+HALFNTH) ,G(N+HALFNTH) , 

&  DIRECDB (N+ HALFNTH) , DIRECRP (N+HALFNTH) 

604  FORMAT ( 3X , F4 . 0 , 9X , ' 0 . 00000 7X - INFINITY 4X , ' -INFINITY ' , 

&  14X,F4.0,2F16.5, F13.5  ) 

END  IF 
END  DO 

320  CONTINUE 

C ****************************************************************************** 
C****************  ************************************************************** 


CL0SE(UNIT-4 , STATUS- ' KEEP ' ) 

CLOSE (UNIT-10 , STATUS- 'KEEP ' ) 

400  CALL  EXIT 
END 
C 
C 

C**************************************************************************** 
C  SUBROUTINE  GAIN  * 

C**************************************************************************** 

SUBROUTINE  GAIN(CDK,CJ,DGAIN,DK,ETT,HKIN,PRAD,SDK,TH,VJ, 

&  WAVM) 

C 

COMPLEX  CJ(l),ETT(l),VJ(l),CQQ,ETHfETD 
DATA  ETA/376.727/ 

CQQ-CMPLX( .0,2 . /ETA) 

CTH«COS( .  01 74533*TH  ) 

STH-SIN(.0174533*TH) 

ETD-(.0, .0) 

IF(TH.GT. 1. ) 

t  ETD-- 60.*( .0,1. )*(CDR-C0S(DK*CTH) ;/(SDK*STH) 

Z.-HK+DK 
ETH-( .0, .0) 

DO  125  I-l.N 
ANG-Z*CTH 

ETT(I)-ETD*CMPLX(COS(ANG) , SIN(ANG) ) 

VJ(I )-CQQ*ETT(I )*VAVM 
IF(PRAD.GT. .0)ETH-ETH+CJ(I)*ETT(I) 

IF(PRAD.LT. .0>CJ(I)-VJ(I) 

125  Z-Z-t-DX 
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IF(PRAD.GT..0)DGAIN-(CABS(ETH)**2)/(30.*PRAD) 

RETURN 

END 

SUBROUTINE  PRCURR 

£*★**★★*★*★*★★★**★****★★★★★**★****★*★★*★**★★**★**★*★*****★***★******■*********** 
C  OUTPUT  STATEMENTS  FOR  CURRENT  DISTRIBUTIONS  * 

C ****************************************************************************** 
VRITE(10,502) 

502  F0RMAT(///17X, 'CURRENT  DISTRIBUTIONS  ON  MONOPOLE') 

WRIT£(10,506) 

506  FORMAT(/5X,'I',12X,'CJ(I)',15X,'CJ(I)',14X,'CJ(I)') 

VRITE(10,508) 

508  F0RMAT(17X, ' (NORM)' |15X, ' (MAG) ' ,8X, ' (PHASE  IN  DEGREES)' ,5X) 
C****************************************************************************** 
c****************************************************************************** 
RETURN 
END 
C 
C 

SUBROUTINE  PRRES(RRAD,RI,XI) 

C 

C  OUTPUT  STATEMENTS  FOR  RADIATION  RESISTANCE,  INPUT  IMPEDANCE,  AND  * 

C  RADIATION  EFFICIENCY  * 

C****************************************************************************** 

C 

REAL*8  ETA 

C 

WRITE( 10, 954)RRAD 

954  FORMAT (IX, 'RADIATION  RESISTANCE  IN  OHMS,  Rrad', 

&  '  (by  integration  of  radiation  pattern)', 

A  5X, 'Rrad  -  '.F10.4  ) 

WRITE( 10 , 952 )RI ,  XI 

952  FORMAT( IX, 'INPUT  IMPEDANCE  IN  OHMS,  Rin  +  jXin' , 

A  42x, 'Rin  -  '.F11.4, 

A  /77X, '  Xin  -  \F11.4  ) 

ETA-RRAD/RI 

WRITE(10,956)ETA 

956  P0RMAT(1X, 'RADIATION  EFFICIENCY,  ETA  -  Rrad/Rin  ', 

A  40X, 'ETA  .  ' , F10.3) 

C****************************************************************************** 

c***********^*****************************************#*********************** 

RETURN 

END 

C 

C 

C 

SUBROUTINE  EARTHTYPE 

C****************************************************************************** 

C  PRINTS  OUT  THE  PROGRAM  NAME  AND  EARTH  TYPE  * 

C****************************************************************************** 

C 

C 

INTEGER  CASE 
CHARACTER*40  TTPEEARTH 

TTPBEARTH  -  '  PERFECT  GROUND' 

CASE  -  1 
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1  WRITE(10, 100) CASE, TYPEEARTB 

100  F0RMAT(///21X,  'PROGRAM  RICHMD6'  ,5X,  'CASE  NO.  '  ,12, ' , ' , (A)) 

C******************************************************^*********************** 

RETURN 

END 

C 

C 

C 

SUBROUTINE  HEADING2(AL,HL,BAR,CMM,IFLAG) 

C  PRINTS  THE  HEADINGS  FOR  OVER  FLAT  LOSSY  EARTH 

C 

C 

C****************************************************************************** 
C  PROGRAM  DESCRIPTION  AND  ECHOED  INPUT  (VITH  EARTH)  * 

C********* ********************* ************************************************ 

c 

CALL  EARTHTYPE 
VRITE(10,224)AL 

224  FORMAT(/ IX, 'MONOPOLE  ELEMENT  RADIUS  IN  WAVELENGTHS,  AL',35X, 

A  'AL  -  ' ,F18. 10) 

WRITE(10,226)HL 

226  FORMAT (IX, 'MONOPOLE  ELEMENT  LENGTH  IN  WAVELENGTHS,  HL',35X, 

A  'HL  .  ' .F18.10) 

WRITE(10,228) 

228  FORMAT (IX, 'RATIO  OF  OUTER  TO  INNER  CONDUCTOR  RADII  OF  THE', 

A  '  COAXIAL' ) 

WRITE(10,230)BAR 

230  F0RMAT(1X, '  LINE  FEED,  BAR' ,62X, 'BAR  »  '.F10.3) 

WRITE(10,83) 

83  FORMAT( IX, 'CONDUCTIVITY  OF  MONOPOLE  ELEMENT  (MEGAMHOS/METER)', 

A  ' ,  CMM  ' ) 

WRITE(10,84)CMM 

84  FORMAT (IX, '  .  -1  FOR  PERFECTLY  CONDUCTING  MONOPOLE  ELEMENT', 

A  30X, 'CMM  «  ' , F10.3) 

WRITE(10,36) 

36  F0RMAT(1X, 'FLAG  FOR  MONOPOLE  ELEMENT  CURRENT  DISTRIBUTIONS', 

A  ',  IFLAG' ) 

WRITE(10,37) 

37  F0RMAT(1X, '  -  0,  MONOPOLE  ELEMENT  CURRENTS  COMPUTED  BY  METHOD', 

A  'OF') 

WRITE(10,11) 

11  FORMAT (7X, 'MOMENTS') 

WRITE(10,38) 

38  F0RMAT(1X, '  -  -1,  MONOPOLE  ELEMENT  CURRENTS  WITH  IDEALIZED', 

A  '  SINUSOIDAL' ) 

WRITE(10, 12) 

12  F0RMAT(8X, 'DISTRIBUTION  AND  FOR  INIFINITE  CONDUCTIVITY  OF  THE' ) 
WRITE(10,13)IFLAG 

13  FORMAT ( 8X , ' MONOPOLE  ELEMENT' ,54X, 'IFUG  -  ‘,14) 

IF(CMM.EQ.-1)GOTO  88 

VRITE(10,86)FMC 

86  FORMAT ( IX , ' FREQUENCY  IN  MEGAHERTZ,  FMC' ,41X, 'PMC  -  ',F8.3) 

88  CONTINUE 

C 

C ****************************************************************************** 

RETURN 


END 


C************ **************************************************************** 

C  SUBROUTINE  CISI  * 

C**************************************************************************** 
SUBROUTINE  Cl SI (Cl, CIN, SI, X) 

O*************************************************************************** 
C  Standard  IBH  Fortran  Subroutine  with  slight  modifications.  * 

C  COSINE  INTEGRAL  AND  SINE  INTEGRAL.  * 

C  X  -  ARGUMENT  (REAL  AND  POSITIVE).  * 

C  Cl  -  Ci(x).  * 

C  SI  -  Si(x) .  * 

C  CIN  -  Cin(x).  * 

O*************************************************************************** 

DATA  GAM, P2/. 57721566, 1.57079632/ 

A-ABS(X) 

IF(A.GT.4. )G0  TO  10 

IF(A.GT..1)G0  TO  3 

IP(A.GT.O. )G0  TO  2 

CI-.O 

CIN- . 0 

SI-.O 

RETURN 

2  X2-A*A 

SI-X*((.03*X2-1. )*X2/18.+1. ) 

CIN-.25*X2'‘((X2/45.-l.)*X2/24.+l.) 

GO  TO  8 

3  Y-(4.-A)*(4.+A) 

SI-X* ( ( ( ( ( 1 . 753141E-9*Y+1 . 568988E-7 ) *Y+1 . 374168E-5 )*Y+6 . 939889E-4 ) 
C*Y+1. 964882E-2 )*Y+4 . 395509E-1 ) 

CIN-  A*A*( ( < ( ( 1 . 386985E-10*Y+1 . 584996E-8 )*Y 

C+l . 725752E-6)*Y+1 . 185999E-4)*Y+4 . 990920E-3 )*Y+1 . 315308E-1 ) 

8  CI-GAM+ALOG(A)-CIN 
RETURN 

10  SI-SIN(A) 

Y-COS(A) 

Z-4./A 

U-( <((<(( (4 . 048069E-3*Z-2 . 279143B-2)*Z+5. 515070E-2 )*Z-7 . 261642E-2 ) 

C*Z+4 . 987716E-2 )*Z-3 . 332519E-3)*Z-2 . 314617E-2 )*Z-1 . 134958E-5)*Z 
C+6.250011E-2)*Z+2.583989E-10 

V-(((((((((-5.108699E-3*Z+2.819179E-2)*Z-6.537283E-2)*Z 
C+7 . 902034B-2 )*Z-4 . 400416E-2 ) *Z-7 . 945556E-3 ) *Z+2 . 601293E-2 )*Z 
C-3 . 764000E-4 )*Z-3 . 122418E-2 )*Z-6 . 646441E-7 )*Z+2 . 500000E-1 
CI-Z*(SI*V-Y*U) 

SI— Z*(SI*UVY*V)*P2 

IP(X.LT..0)SI— SI 

CIN»GAM+ALOG(A)-CI 

RETURN 

END 


C***************************************************************************** 
C  SUBROUTINE  FRILLS  * 

C*********************************** ****************************************** 

SUBROUTINE  FRILLS(AR,BAR,DK,NEQ,VJ) 

C  FRILLS  acts  up  the  voltage  column  VJ(I).  * 

C  VJ(I)  -  voltage  column  for  perfectly  conducting  vire  dipole  in  * 
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non 


20 


fr«*  space  using  Calerkin's  aethod  and  sinusoidal  bases,  and  * 

aatching  the  boundary  conditions  on  the  surface  of  the  ir  .  * 

c****J*^*^*****^i**"********#‘*"***;‘*********;***************^****** 

^^DKi;«?GI(20),VJ(l),GII,QST,WST 

DATA  PI , TP/3 . 14159265359 , 6 . 28318530718/ 

IDH-20 

DO  20  I»1,NEQ 
VJ(IM.0,.0) 

VJ(l)-{l.,-0) 

IF(BAR.LE.1.)RETIRN 

VJ(1)«(.0,.0) 

NSV-NEQ+1 

SDK-SIN(DK) 

CDK-COS(DK) 

BAL-ALOG(BAR) 

QST-CMPLX ( . 0 , 1 . / ( 4 . *BAL* SDK ) ) 

BK-AK*BAR 
AKS-AK*AK 
BKS-BK*BK 
LIM-NSW+1 

IF(LIM.GT.IDM)LIH-IDM 
NFB-6 

NPB-2*(NPB/2) 

NPP-NPB+1 
PHA-. 0174533*20. 

DPH-PHA/NPH 


PH-.O 

DO  90  LPB-1,2 
WST-DPH*QST/(3.*PI) 

SGH— 1. 

DO  80  IPB-l.NPP 

WF-3.+SGN 

IF(IPB.E0.1)VF.l. 

IF(IPB.EQ.NPP)VF-1. 

CPB-COS(PB) 

IF(IPB.GT.1)G0  TO  40 
IF(LPB.GT. l)GO  TO  40 
CPH-C0S(DPB/10. ) 

40  RSI-2. *AKS*(1.-CPB) 

RS2-AKS+BKS-2 . *AK*BK*CPH 

RHl-DSQRT(RSl) 

RH2«DSQRT<RS2) 

CALL  CISI(CA,CIN,SA,RH1) 
CALL  CISI(CB,CIN,SB,RH2) 
GI(l)-2 .*CMPLX(CB-CA, SA-SB) 
DO  50  1-2, LIM 
DZ-DK*(I-1 ) 

DZS-DZ*DZ 

RA-DSQRT(RS1+DZS) 

RB-DSQRT(RS2+DZS) 

CALL  CISI(C1,CIN,S1,RA*DZ) 
CALL  CISI<C2,C»,S2,RB+DZ) 
CP-OffU(C2-Cl,Sl-S2) 
RAM-RS1/(RA+DZ) 

RBH-RS2 / (RB+DZ) 

CALL  CISI(CI,CIN,S1,RAH) 
CALL  CISI(C2,CIN,S2,RBH) 
GH-OffU(C2-Cl,Sl-S2) 
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BGZ-CMPLX ( COS( DZ ) , SIN ( DZ ) ) 

50  GI(I)-GP*BGZ+GM/BGZ 

V J ( 1 ) -VJ ( 1 ) ♦WF*VST* ( GI ( 2 ) -CDK*GI ( 1 ' ) 

IF(NEQ.LE. 1)G0  TO  78 
Kl-0 

SO  60  I -2, NEQ 
Kl-Kl+1 
K2-KU1 
K3-R2+1 

IF(K3.GT.IDM)G0  TO  60 
GP-G1(K1)-2.*CDK*GI(K2)+GI(K3) 

V J ( I ) -V J ( I ) +VF*VST*GP 
60  CONTINUE 
78  SGN--SGN 
80  PH-PH+DPH 

DPH-(PI-PHA)/NPH 
90  PH-PHA 

VJ(1)-2.*VJ<1) 

RETURN 

END 

C 

C 

C***************************************************************************** 


C  SUBROUTINE  TPLZ  * 

C***************************************************************************** 

SUBROUTINE  TPLZ(C,U,VJ,V,Z,IBR,IVR,I12,N£Q) 

C********************************************************************* ******** 


C  MODIFIED  VERSION  OF  SUBROUTINE  FURNISHED  BT  CHARLES  KLEIN.  * 
C  SOLVES  SIMULTANEOUS  LINEAR  EQUATIONS.  * 
C  SET  ItfR  -(1  OR  0)  TO  GET  (PRINTOUT  OR  NO-PRINTOUT).  * 
C  SET  112  -  1  ON  FIRST  CALL,  WHERE  MATRIX  INVERSION  IS  REQUIRED.  * 
C  SET  112  -  2  IF  MATRIX  Z  HAS  ALREADY  BEEN  INVERTED  ON  PREVIOUS  CALL.  * 
C  NEQ  .  NUMBER  OF  SIMULTANEOUS  UNBAR  EQUATIONS.  * 
C  Z(J)  IS  THE  FIRST  ROW  OF  THE  TOEPUTZ  IMPEDANCE  MATRIX  * 
C  VJ(J)  -  INPUT  VOLTAGE  COLUMN  * 
C  C(J)  -  OUTPUT  CURRENT  COLUMN  * 
C  U(J)  AND  V(J)  ARE  WORK  ARRAYS  OF  LENGTH  NEQ  * 
C  IF  IER  -  0  ,  NO  ERROR  OCCURRED  * 


C***************************************************************************** 


COMPLEX  C(1),U(1),VJ(1),V(1),Z(1) 

COMPLEX  ALMDAIALPHA,C1,C2,C0EF,FAC,TAU1,V,V1,V2 
2  FORMAT(1X,I5,7X,F10.4,7X,F15.7, 7X.F10.1) 

7  FORMAT (5X) 

IF(NBQ.GT.1)G0  TO  8 
C(1).VJ(1)/Z<1) 

CN0R-CABS(C(1)) 

GO  TO  100 

8  IF(I12.NR.1)G0  TO  45 
N-NEQ-1 

IBR-0 


C  NORMALIZE  INPUT  MATRIX 


* 


TAUl-Z(l) 

DO  10  IX-1,N 
10  Z(II)-Z(1I>1)/TAU1 
ALMDA-1.-Z(1)*Z(1) 
U(1)-Z<1) 

1-2 

15  KK-I-1 


C-ll 


ALPHA- (.0,.0) 

DO  20  H-1'KK 
LL-I-M 

20  ALPHA-ALPHA+U(M)*Z(LL) 

ALPHA—  (ALPHA«Z(I)) 

IF(CABS( ALPHA). EQ. . 0)G0  TO  130 

COEP-ALPHA/ALMDA 

ALMDA-ALMDA-COEF*ALPHA 

DO  30  J-1,KK 

UI-J 

30  V(J)«U(J)+COEF*U(L) 

DO  40  J-l.KK 

40  U(J)-V(J) 

U(I)*COEP 
IF(I .GE.N)GO  TO  45 
1-1*1 
GO  TO  15 

£****************************#************************************************ 

C  THE  FOLLOWING  COMPUTES  THE  ELEMENTS  OP  THE  INVERSE  * 

£1^**^****%**************************##*****#********************************* 

45  NH-(NE0+l)/2 

FAC-ALMDA*TAU1 
NP-NEQ+1 
CNOR-.O 
DO  90  I-l.NH 
IF(I.NE.1)G0  TO  55 
W(1)-1./FAC 
DO  50  J-2,NEQ 

50  V(J)«U(J-1)/FAC 
GO  TO  70 

55  Cl-U(I-l) 

NPI-NP-I 

C2-U(NPI) 

DO  60  JJ-l.N 

J-NP-JJ 

NPJ-NP-J 

60  V(J)-V(J-1)*<C1*U(J-1)-C2*U(NPJ))/FAC 
W(1)-U(I-1)/PAC 

C  MATRIX  MULTIPLY  * 

C***************************************************************************** 

70  V-(.0,.0) 

Vl-(.0,.0) 

DO  80  J-1,NEQ 
V2-VJ(J) 

V-V*V2*W(J) 

80  V1-VUV2*W(NPJ) 

C(I)-V 

NPI-NP-I 

C(NPI)-V1 

XF(IVR. LE.O)GO  TO  90 
CA-CABS(V) 

IF(CA.GT.CNOR)CNOR-CA 

CA-CABS(Vl) 

IF(CA.GT.CNOR)CNOR-CA 

90  CONTINUE 

100  IP(IWR.LE.O)GO  TO  120 


CATMINT* OUT  THE  SOLUTION  FOR  THE  CURRENTS  C(J)  * 
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£★★★★★★*★★★*★*★★★★*★**★★★★★★★★*★**★★*★★★***★**★*★★*★***★***★**★*************** 
IF(CNOR.LE.O. )CN0R.l. 

DO  I-l,(NEQ+l)/2 
V-C(I) 

CA-CABS(V) 

CN-CA/CNOR 

PH-.O 

IF(CA.GT .0. )PH«57 . 29578*ATAN2(AIMAG(V) , REAL(V) ) 

IF(I.NE. ((NEQ+1)/2))VRITE(10,2)I,CN, (CA*2),PH 
IF(I.EQ. ((NEQ+1)/2))VRITE(10,22)I,CN, (CA*2) ,PH 
22  F0RMAT(1X,I5,'  (base)' ,fl0.4,7x,fl5.7,7x,fl0.1) 

END  DO 
120  RETURN 
130  IBR-I 
RETURN 
END 
C 
C 


C  SUBROUTINE  TSPAR  * 

C***************************************************************************** 
SUBROUTINE  TSPAR(AK,DK,NEQ,Z) 

C**************** ************************************************************* 
C  TSPAR  sets  up  the  impedance  matrix  Z(J).  * 

C  Z(J)  -  First  rov  of  impedance  matrix  for  perfectly  conducting  * 

C  thin-vlre  dipole  in  free  space,  using  Galerkin’s  method  vith  * 

C  overlapping  sinusoidal  basis  functions  and  matching  the  * 

C  boundary  conditions  on  the  surface  of  the  vire.  * 

C  AK  -  k*a,  vhere  k  «  2*pi /lambda  and  a  •  vire  radius.  * 

C  DK  -  k*d,  vhere  d  -  segment  length.  * 

C  NEQ  -  number  of  simultaneous  linear  equations.  * 

C***************************************************************************** 

REAL*8  DZS.RS 

COMPLEX  EID(90) ,EM(90) ,EP(90) ,Z(1) 

COMPLEX  CEM , CEP , EMD , EPD , EMD2 , EPD2 ,Z11,Z22,G11,Q11 
DIMENSION  CID(90) , SID(90) ,CM(90) ,CP(90) , SM(90) , SP(90) 

DATA  GAM.P2/. 577215664,1. 57079632/ 

DATA  ETA.PI/376.727.3.14159/ 

IDM-90 

1  FORMAT ( 3X , ' MUST  INCREASE  DIMENSIONS  IN  SUBROUTINE  TSPAR') 

2  FORMAT (3X, 'ACTUAL  DIMENSION  I DM  -  ',I5,6X, 

2 'REQUIRED  DIMENSION  MAX2  -  ',15) 

IF(NEQ.LE.O)R£TURN 
MAX2-NEO+2 
DO  14  I.l.NEQ 
14  Z(I)-(.0,.0) 

IF(MAX2.LE.IDH)G0  TO  16 
WRITE(lO.l) 

VRITE(10,2)ZDM,MAX2 

RETURN 

16  TDK-2. *DK 
S11-.0 
S13-TDK 
S21.DK 
S23-3.*DK 
DO  20  N.1.MAX2 
I-N-l 
DZ-I*DK 
CID(N)-COS(DZ) 

SID(N)-SIN(DZ) 
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20  EID(N)-CMPLX(CID(N) , SID(N) ) 

CDK-COS(DK) 

SDK-SIN(DK) 

EPD-CHPLX(CDK, SDK) 

EMD-CMPLX(CDK, -SDK) 

EPD2-EPD*EPD 
EHD2  «EMD*EMD 
CEM-2 . *CDK+EMD 
CEP-2 . *CDK+BPD 
AK2-AK*AK 

CSS-ETA/ (8 . *PI *SDK*SDK) 

NPB-6 

NPB-2*(NPH/2) 

NPP-HPB+1 
PHA-. 0174533*20. 

DPB-PHA/NPH 

PB-.O 

DO  100  JPB-1,2 

CST-DPB*ETA/ ( 24 . *PI*PI*SDK*SDK) 

C22-DPB/<3.*PI) 

SGN— 1. 

DO  80  ZPH-1,MPP 
CPB-COS(PB) 

SPH-SIN(PB) 

IF(IPH.GT. 1 )G0  TO  30 
IF( JPB.GT . 1)G0  TO  30 
PB0-DPB/10. 

CPH-COS(PBO) 

SPH-SIN(PBO) 

30  RH-AK*SPH 

RS-2 . *AK2*( 1 . -CPB) 

RK-DSQRT(RS) 

WF-3.+SGN 

IF(IPB.EQ.1)VF-1. 

1F(IPB.E0.NPP)WF-1. 

UST-WF*CST 

V22-VF*C22 

DO  40  N-l ,MAX2 

I-N-l 

DZ-I*DK 

DZS-D2*DZ 

R-DSQRT ( RS+DZS ) 

ARG-R+DZ 

IF(N.EQ.  DARC-RK 

CALL  CISI(CP(N),CIN,SP(N),ARG) 

EP(N)-CMPLZ(CP(N) , -SP(M) ) 

IF(N.GT. 1)G0  TO  38 
CM(1)-CP<1) 

SM(1)-SP(1) 

EM(1)-EP(1) 

GO  TO  40 
38  ARG-RS/ARG 

CALL  CISI(CM(N)tCINtSM(N),ARG) 

EM(H)-CNPLZ(CM(N) , -SM(N) ) 

40  CONTIHUB 

R-4.*(-CM(2)+2.*CP(l)-CP(2)) 

A+2.*C1D(3)*(4CR(3)-2.*CM(2)42.*CP(1)-2.*CP(2)*CP(3)) 

»42.*SID(3)*(-S«(3)42.*SII(2)-2.*SP(2)4SP(3)) 

X-4.*(SM(2)-2.*SP(1)4SP(2)) 

C+2.*CTD(3)*(-SR(3)42.*SI!C2)-2.*SP(1)42.*SP(2)-SP(3)) 


o  o 


D+2.*SID(3)*(-CN(3)+2.*CH(2)-2.*CP(2)+CP(3)) 

Z(1)»Z(1)+VST*CMPLX(R,X) 

IP(N£Q.EQ. 1)G0  TO  70 

R»2.*CID(2)*(-CM(3)+3.*CH(2)-4.*CP(1)+3.*CP(2)-CP(3)) 

B*2.*SID(2)*(+SN(3)-2.*SM(2)+2.*SP(2)-SP(3)) 

F+CID(4)*(+CM(4)-2.*CM(3)+CM(2)+CP(2)-2.*CP(3)+CP(4)) 

G+SID(4)*(-SH(4)+2.*SM(3)-SM(2)+SP(2)-2.*SP(3)+SP(4)) 

X-2.*CID(2)*(SM(3)-3.*SM(2)+4.*SP(1)-3.*SP(2)+SP(3)) 

H  ♦2.*SID(2)*(CM(3)-2.*CM{2)+2.*CP(2)-CP(3)) 
I-fCID(4)*(-SH(4)+2.*SM(3)-SM(2)-SP(2)+2.*SP(3)-SP(4)) 
J+SID(4)*(-CM(4)+2.*CM(3)-CM(2)+CP(2)-2.*CP(3)+CP(4)) 
Z(2)-Z(2)+tfST*CMPLX(R,X) 

IF(NEQ.EQ.2)G0  TO  70 
Sl-DK 

DO  60  N-3.NBQ 

Ml-N-1 

M2-N-2 

Nl-N+1 

N2-N+2 

CPA-CP(M2)-2.*CP(M1)+CP(N) 

CPB-2.*CP(N)-CP(M1)-CP(N1) 

CPC-CP(N2)-2.*CP(N1)*CP(N) 

CMA-CM(M2)-2.*CM(M1)+CM<N) 

CMB-2 . *CM(N)-CM(N1 )-CM(Ml ) 

CMC-CM(N2 )-2 . *CM(N1 )+CM(N) 

SPA-SP(M2)-2.*SP(M1)+SP(N) 

SPB-2.*SP(N)-SP(M1)-SP(N1) 

SPC-SP(N2 )-2 . *SP(N1 )+SP(N) 

SMA«SM(H2)-2.*SM(M1)+SM(N) 

SMB-2. *SM(N)-SM<N1)-SM(M1) 

SMC-SM(N2)-2 .*SM(N1)+SM(N) 

R-CID(M2)*(CPA*CMA)+2.*CID(N)*(CPB*CMB)+2.*SID(N)*(SPB-SMB) 
K  +CID( N2 ) * ( CPC+CMC )+SID(N2)*( SPC- SMC ) 
IF(N.GT.3)R»R+SID(M2)*(SPA-SMA) 

X— CID(M2)*(SPA+SMA)-2.*CID(N)*(SPB+SMB)+2.*SID(N)*(CPB-CMB) 
L  -CID(N2)*(SPC+SMC)+SID(N2)*(CPC-CMC) 
IF(N.GT.3)X-X+SID(M2)*(CPA-CMA) 

60  Z(N)-Z(N)+VST*CMPLX(R,X) 

70  PB-PB+DPB 
80  SGN— SGN 

DPB-(PI-PBA)/NPB 
100  PB-PBA 
RETURN 
END 


C******************** ********************************************************* 

C  SUBROUTINE  ZSURF  * 

C************************"***************************************************** 


SUBROUTINE  ZSURF (AX, CMM, PMC, ZS) 


C  ZSURF  calculates  the  surface  lnpedance  ZS  for  thin  vire. 

C  AX  -  k*a  vhere  k  •  2*pi/laabda  and  a  -  vire  radius. 

C  CMM  -  conductivity  of  vire  (negaahos/neter) 

C  CMM  -  -1  for  perfect  conductivity. 

C  FNC  -  frequency  (eegaherts). 


COMPLEX  BBS,BES1,ZS 

DATA  ETA, SQT , TP/376 . 72727 , 1 . 41421356 , 6 . 2831853/ 
SQSVE.l . E6*SQRT(CHM /TP/PMC/8. 85433) 
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X- AR*SQSVE 
IF(X.GT.8. )G0  TO  50 
T-X/8. 

T2«T*T 

T4«T2*T2 

BER«(( ((<(-. 901£-5*T4+ . 122552E-2 >*T4- . 08349609 )*T4 
2+2. 641914)*T4-32. 363456 )*T4+113. 77778 )*T4-64. )*T4+1. 
BEI»( (((((. 11346E-3*T4- . 01103667>*T4+ . 52185615)*T4 
2-10. 567658)*T4+72. 817777 )*T4-113.77778)*T4+16. )*T2 
BERP-X*T2* ((<(((-. 394E-5*T4+ . 45957E-3 ) *T4- . 02609253 )*T4 
2+ . 66047849 )*T4-6 . 0681481 )*T4+14 . 222222)*T4-4 . ) 

BEIP.X*( <((((. 4609E-4*T4- . 379386E-2 )*T4+ . 14677204)*T4 
2-2. 3116751 )*T4+11. 377778 )*T4-10. 666667 )*T4+. 5) 

BES-CMPLX ( BER , BEI ) 

BES1. . 707107*CHPLX(BERP-BEIP, BERP+BEIP) 

GO  TO  100 
50  XP-.70710681*X 

XI- l./X 

F-( (- .0459205*X1+ . 390625E-2 )*X1+ .08838835)*X1+1 . 

T«( (- . 04603559*X1- . 0625)*X1- . 08838835)*X1- . 39269907+XP 
BES-P*CMPLX(COS(T) ,SIN(T>) 
F-((.11290231*X1+.03515625)*X1-.26516505)*X1+1. 

T-( ( . 1160097*X1+ . 1875)*X1+ . 26516505 )*X1+1 . 1780972+XP 
BES1-F*CKPLX(C0S(T) , SIN(T) ) 

100  ZS— CMPLX( 1 . , -1 . )*ETA*BES/BES1/SQT/SQSWE 
RETURN 
END 
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APPENDIX  D 

COMPUTER  PROGRAM  WAIT-SURTEES  FOR  THE  INPUT  IMPEDANCE  OF  A 
MONOPOLE  ELEMENT  ON  A  DISK  GROUND  PLANE  ABOVE  FLAT  EARTH 
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COMPUTER  PROGRAM  WAIT-SURTEES 
by 

JACK  11.  RICHMOND 
December  29,  1989 


INTRODUCTION  1 

Appendix  1  presents  Richmond's  computer  program  WAIT-SURTEES. FOR 
together  with  the  subroutines  C1SI,  FRILLS,  TPLZ,  TSPAR  and  ZSURF. 
This  FORTRAN  program  calculates  (lie  impedance  of  a  vertical  monopole 
antenna  centered  on  a  circular  discover  the  Hat  lossy  earth.  This  program 
combines  the  following: 

a)  Richmond’s  moment  method  for  the  impedance  Zx  of  a  vertical 
monopole  antenna  on  an  infinite  ground  plane* and 

b)  The  theory  of  Wait  and  Surtees  for  the  change  A Z  in  the  antenna 
impedance,  where  A Z  =  Zt  -  Z ^  and  Z/  denotes  the  impedance  of  the 
vertical  monopole  on  a  finite  circular  ground  plane  over  the  flat  earth. 

See:  (J.  R.  Wait  and  W.  J.  Surtees,  "Impedance  of  Top- Loaded  Antenna 
of  Arbitrary  Length  Over  a  Circular  Grounded  Screen,"  J.  Appl.  Phvs., 
Vol.  25,  pp.  553-555,  May  1954). 

Comment  statements  have  been  inserted  in  the  main  computer  program 
and  in  each  subroutine  to  assist  the  user.  Only  a  few  brief  additional 
comments  will  be  required  in  this  Introduction. 

In  calculating  Zw,  the  monopole  is  divided  into  segments  of  equal  length 
and  the  unknown  current  distribution  is  expanded  in  overlapping  sinu¬ 
soidal  basis  functions.  Thus,  l(z)  is  taken  to  be  piecewise  sinusoidal.  The 
magnetic-frill  model  is  employed  (rather  than  the  slice-generator  model), 
and  boundary  matching  is  enforced  on  the  surface  of  the  monopole  rather 
than  on  the  axis.  The  wire  radius  is  assumed  to  be  much  smaller  than 
the  wavelength.  The  wire  monopole  may  be  assigned  perfect  conductiv¬ 
ity  or  finite  conductivity  as  desired.  With  Galerkin’s  method,  the  calcu¬ 
lated  impedance  Z,*,  is  believed  to  be  accurate  for  short,  medium  and  long 

1  Appreciation  is  expressed  to  The  MITRF.  Corporation  for  sponsoring  this  report  . 

The  compoler  program  WAIT-SURTEES. FOR  was  developed 
by  Richmond  in  1979  with  other  sponsorship. 

*  of  Infinite  conductivity 
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monopoles.  If  the  monopole  length  exceeds  6  wavelengths,  however,  one 
may  wish  to  increase  the  dimensions  (1DM=99)  in  the  main  program  and 
subroutine  TSPAR. 

In  the  theory  of  Wait  and  Surtees,  the  mor.opole  is  assumed  to  have  a 
sinusoidal  current  distribution.  (It  <1  oes  »/>(  appear  difiirull  to  generalize 
this  to  a  piecewise-sinusoidal  distribution,  but  we  have  not  attempted  this.) 
Since  the  current  distribution  departs  significantly  from  the  sinusoidal  lorni 
when  the  monopole  length  exceeds  one-half  wavelength,  A Z  and  Zj  may 
begin  to  lose  reliability  as  the  monopole  length  increases.  We  have  not 
investigated  this  possible  problem. 

Appendix  II  presents  the  output  data  generated  by  WAIT-SURTEES.FOR 
on  a  VAX  computer  with  the  same  input  data  indicated  in  Appendix  1.  This 
output  shows  excellent  agreement  with  the  original  results  obtained  in  1979 
on  a  DATACRAFT  computer.  This  indicates  that  no  additional  double- 
precision  operations  are  required  for  VAX  operation. 

Richmond  has  shown  that  WAIT- SURTEES. FOR  is  useful  even  for 
a  monopole  antenna  on  a  circular  disk  in  free  space.  [J.  II.  Richmond, 
“Monopole  Antenna  on  Circular  Disk,”  IEEE  Trans.,  Vol.  AP-32,  pp.  1282- 
1287,  December  1984).  For  this  case,  set  ER=1  and  SIC=0  in  the  main 
program. 

In  WAIT-SURTEES.FOR  the  monopole  is  centered  on  a  circular  disk 
which  may  lie  on  the  surface  of  the  earth,  or  it  may  be  located  any  distance 
above  the  earth  surface.  In  the  main  program  1IDL  denotes  the  height  of 
the  circular  disk  above  the  flat  earth,  measured  in  free-space  wavelengths. 

For  a  monopole  on  a  circular  disk  on  the  surface  of  the  earth,  Rich¬ 
mond  has  shown  satisfactory  agreement  between  WAIT-SURTEES.FOR 
and  Richmond’s  moment  method  (which  enforces  the  boundary  conditions 
to  determine  the  current  distributions  on  the  monopole  and  the  disk).  See: 

(J.  H.  Richmond,  “Monopole  Antenna  on  Circular  Disk  Over  Flat  Earth,” 
IEEE  Trans.,  Vol.  AP-33,  pp.  033-637,  June  1985.] 

In  the  output  data  of  Appendix  II,  the  resistance  RFIN  is  plotted  as  the 
dashed-line  curve  of  Figure  5  in  (Richmond,  1985].  The  reactance  XF1N  in 
Appendix  II  should  have  been  plotted  as  the  dashed-iine  curve  of  Figure  0 
in  (Richmond,  1985].  By  mistake,  however,  the  dashed-iine  curve  of  Figure 
6  shows  the  output  of  WAIT-SURTEES.FOR  for  a  monopole  on  a  circular 
disk  in  free  space. 
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Appendix  I.  WA IT-SURT5ES .FOR 


WAIT-SURTEES.  FOR  WMT.l 

impedance  of  monopole  at  center  of  circular  disk  on  flat  earth. 

SEE:  WAIT  AND  SURTEES,  "IMPEDANCE  OF  TOP-LOADED  ANTENNA  OF 

ARBITRARY  LENGTH  OVER  A  CIRCULAR  GROUNDED  SCREEN,"  J.  APPL.  PHYS . , 
VOL.  25,  PP.  553-555,  MAY  1954. 

LINK  Cl  SI ; TRILLS; TPLZ ; TSPAR; ZSURF 

COMPLEX  COO, CJI , CJA, CJB, DELZ 

COMPLEX  EC, EGZ, EJH, EJS, ETA2, ETD, ETH 

COMPLEX  GM,GP,Q1,02,C3,RC 

COMPLEX  Y1 1 , Zll , ZFIN, ZINF , ZS, ZSG 

COMPLEX  C J  (99)  ,  GI  (99)  ,U(99)  ,  VJ  (99)  ,  W  (99)  ,  ZJ  (99) 

DATA  ETA, PI, TP/376. 730366239, 3. 14159265359, 6.28318530718/ 

DATA  EO, UO/B . 85418533677E-12, 1 .25663706144E-6/ 

DATA  P2/1 .57079632679/ 

DATA  IDM/99/ 

AI ,  AL,  AM  -  RADIUS  OF  MONOPOLE  WIRE  IN  ( INCHES ,  WAVELENGTHS ,  METERS)  . 
BAR  -  RATIO  OF  OUTER  RADIUS  AND  INNER  RADIUS  OF  COAXIAL  FEED. 

BL.BM  -  OUTER  RADIUS  OF  CIRCULAR  DISK  (WAVELENGTHS,  METERS)  . 

CMM  -  CONDUCTIVITY  (MEGAMHOS/METER)  OF  MONOPOLE  WIRE. 

CMM  -  -1.  FOR  PERFECTLY  CONDUCTING  MONOPOLE. 

ER  -  RELATIVE  PERMITTIVITY  OF  EARTH. 

FMC  -  FREQUENCY  (MEGAHERTZ)  . 

HDL  -  HEIGHT  OF  CIRCULAR  DISK  ABOVE  THE  FLAT  EARTH  (WAVELENGTHS)  . 

HDL  -  -1.  TOR  CIRCULAR  DISK  ON  THE  SURFACE  OF  THE  EARTH. 

HL,HM  -  LENGTH  OF  MONOPOLE  (WAVELENGTHS ,  METERS) 

SIG  -  CONDUCTIVITY  OF  EARTH  (MHOS/METER) . 

WAVM  -  WAVELENGTH  IN  FREE-SPACE  (METERS)  . 

2  FORMAT (IX, 8F11 . 5) 

5  FORMAT (1H0) 

AL* .003 
BAR  -  3. 

CMM— 1  . 

ER-4  . 

FMC-300. 

HL* .229 
SIG*. 001 
WAVM-300 . /FMC 
IWR-l 

OMEG-TP*FMC«l ,E6 
EC-CMPLX  (ER,  -SIG/  (OMEG*E0)  ) 

ETA2-ETA/CSQRT (EC) 

AK«TP"AL 
HK-TP*HL 
SK«2  .*HK 

EJS-CMPLX  (COS  (SK)  ,  SIN  (SK)  ) 

CHK-COS (HK) 

SHK-SIN(HK) 

EJH -CMP LX  (CHK,  SHK) 

NS-15.*HL 

IF (NS. LT. 6  INS-6 

IF (NS . GT . IDM) NS-IDM 

IG-NS/2 

NS-2*IG 

N-NS-1 

ZS»( .0, .0) 

IF (CMM . GT . 0 . ) CALL  ZSURF (AK, CMM, FMC, ZS) 

ZS  -  SURFACE  IMPEDANCE  OF  MONOPOLE  WIRE. 

DK-HX/IG 
CDK-COS (DK) 

SDK-SIN (DK) 

ZJ(N)  -  FIRST  ROW  OF  IMPEDANCE  MATRIX  FOR  DIPOLE  IN  TREE  SPACE, 

USING  GALERXIN'S  METHOD  WITH  OVERLAPPING  SINUSOIDAL  BASIS  FUNCTIONS. 
(DIPOLE  LENGTH  -  TWICE  THE  MONOPOLE  LENGTH,  USING  IMAGE  THEORY.) 

CALL  TSPAR  (AK,  DK,N,  ZJ) 

IF (CMM . LT . 0 . ) GO  TO  52 
GK-2 .* (PK-CDK"SDK) 
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GL-SDK-DX«CDK 
FH-4 . *PI*AX*SDK*SDK 
ZJ (1 ) -ZJ  < 1 ) +ZS*GX/FH 
ZJ(2)-ZJ(2>+ZS*GL/FH 
52  112-1 

C  SET  UP  THE  MOMENT-METHOD  VOLTAGE  COLUMN  VO  (M)  FOR  DIPOLE  IN  FREE  SPACE. 

CALL  FRILLS  (AX,  BAR,  DK.N.CJ) 

DO  60  I— 1 , N 
K— 1+IABS (IG-I) 

60  VJ(I)-CJ(K) 

C  SOLVE  THE  SIMULTANEOUS  LINEAR  EQUATIONS  TO  DETERMINE  THE  CURRENTS  CJ(N) 

C  ON  THE  WIRE  DIPOLE  IN  FREE  SPACE.  (THE  IMPEDANCE  MATRIX  IS  TOEPL1TZ.) 

CALL  TPLZ (CJ, U, VJ , W, ZJ , IER, INR, 112, N) 

C  ZINF  -  IMPEDANCE  OF  MONOPOLE  ANTENNA  OVER  INFINITE  GROUND  I’LANE . 

ZINF- . 5/CJ (1G) 

WRITE ( 6 , 2 ) AL, HL, ZINF 
WRITE  (15, 2)  AL,  HL,  ZINF 
WRITE  (6, 5) 

WRITE (15, 5) 

HDL— 1 . 

HDX-TP*HDL 

ZSG-ETA2 

IF(HDL.LE. .0)GO  TO  70 
RC-  (ETA2-ETA)  /  (ETA2+ETA) 

EJH-CMPLX (COS (HDK) , SIN (HDX) ) 

ZSG-ETA* (EJH+RC/E JH) / (EJH-RC/EJH) 

70  CONTINUE 

DO  100  IBL-1,16 

BL«.1*IBL 

BK«TP*BL 

RK-SQRT (BK*BK+HX*HK) 

C  ZINF  -  IMPEDANCE  OF  MONOPOLE  ON  INFINITE  CIRCULAR  DISK. 

C  ZFIN  -  IMPEDANCE  OF  MONOPOLE  AT  CENTER  OF  FINITE  CIRCULAR  DISK. 

C  DELZ  -  ZFIN  -  ZINF. 

C  CALCULATE  DELZ  USING  THE  FORMULA  OF  WAIT  AND  SURTEES. 

CALL  CIS1 (CP , CIN, SP,  2 . * (RX+HK) ) 

CALL  CISI(CM,CIN,SM,2.*(RK-HK)) 

Ql-CMPLX (CP, P2-SP) *EJS+CKPLX (CM, P2-SM) /EJS 
CALL  Cl SI (CP,CIN,SP,RK+BK+HK) 

CALL  CISI (CM,CIN,SM,RK+BK-HX) 

CALL  CISI (CB, CIN, SB, RX+BX) 

Q2-CMPLX (CP , P2-SP ) *E JH+CMPLX (CM, P2-SM) /EJH-CMPLX (CB, P2-SB) 

CALL  CISI (Cl, CIN, SI, 2 .*BK1 
Q3-CMPLX (Cl , P2-SI ) 

DELZ-Q1-4 . »CHK*Q2+2 . *CHX*CHX*Q3 
DELZ— ZSG* DELZ/ (4 . «PI*SHK*SHX) 

ZFIN-ZINF+DELZ 

WRITE (15, 2) BL, ZFIN 
100  WRITE (6, 2) BL, ZFIN 
400  CALL  EXIT 
END 
C 
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SUBROUTINE  CISI (CI,CIN, SI,X) 

Standard  IBM  Fortran  Subroutine  with  alight  modifications. 

COSINE  INTEGRAL  AND  SINE  INTEGRAL. 

X  «  ARGUMENT  (REAL  AND  POSITIVE) . 

Cl  -  Ci (X) . 

SI  -  Si (X) . 

CIN  -  Cin(x) . 

DATA  GAM, P2/. 57721566, 1 .57079632/ 

A-ABS (X) 

IF (A. GT. 4.) GO  TO  10 

IF (A . GT . . 1) GO  TO  3 

IF (A .GT . 0 . ) GO  TO  2 

CI-.  0 

CIN-.O 

SI-.O 

RETURN 

2  X2— A*A 

SI-X« ( <.03*X2-1 .) -X2/18.+1 . ) 

CIN-.25«X2*((X2/«5.-l.)*X2/24.+l.) 

GO  TO  8 

3  Y-(4.-A)*(4.+A) 

SI-X* ( ( ( ( <1 ,7S3141E-9*Y+1 .56B9BSE-7) *Y+1 . 374168E-5) *Y+6 . 939889E-4) 
C*Y+1 .964882E-2) *Y+4 . 395509E-1 ) 

CIN-  A*A« ( ( ( ( (1 .3B6985E-10*Y+1 . 584996E-B) *Y 

C+1.725752E-6)*Y+1 .185999E-4) *Y+4 . 990920E-3) *Y+1 . 315308E-1) 

B  CI-GAM+ALOG (A) -CIN 
RETURN 

10  SI-SIN  (A) 

Y-COS (A) 

Z— 4 . /A 

U- ((((((( (4 . 048069E-3*Z-2 .279143E-2) *Z+5 . 515070E-2) *Z-7 .261642E-2) 
C*Z+4 . 98771 6E-2) *Z-3 . 332519E-3) *Z-2 . 314617E-2) *Z-1 . 134958E-5)  »Z 
C+6. 25001 1E-2)*Z+2.5839B9E-10 
V- (((((((( (-5 . 10B699E-3»Z+2 . 819179E-2) «Z-6.537283E-2) *Z 
C+7 . 902034E-2) *Z-4 . 400416E-2) »Z-7 . 945556E-3) »Z+2 . 601293E-2) *Z 
C-3 . 764000E-4) »Z-3.122418E-2) «Z-6 . 646441E-7) *Z+2 . 500000E-1 
CI-Z* (SI*V-Y*U) 

SI  — Z*  (SI*U+Y*V)  +P2 
IF  (X .  LT .  ,0)SI  — SI 
CIN-GAM+ALOG  (A)  -CI 
RETURN 
END 
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FRILLS.] 


SUBROUTINE  FRILLS (AK,  BAR, DK, NEQ, VJ) 
c  FRILLS  sets  up  the  voltage  column  VJ ( I ) . 

C  VJ(I)  -  voltage  column  for  perfectly  conducting  wire  dipole  in 
C  free  apace  uaing  Galarkin'a  method  and  ainusoidal  bases,  and 

C  matching  the  boundary  conditions  on  the  surface  of  the  wire. 

C  Using  magnetic-frill  model  for  center-fed  dipole. 

REAL* S  DZS , RSI , RS2 

COMPLEX  EGZ, GM, GP, GI (20) ,  VJ (1 ) , GI I , QST, WST 
DATA  PI,  TP/3. 1415 9265359,  6. 2831 853011 8/ 

IDM-20 

DO  20  I-l.NEQ 
20  VJU)-(.0,.0) 

VJ<1)-(1.,  .0) 

IF (BAR ■ LE . 1 . ) RETURN 
VJ(l)-(.0,  .0) 

NSW-NEQ+ 1 
SDX-SIN(DX) 

CDK-COS (DK) 

BAL-ALOG  (BAR) 

QST-CMPLX ( . 0, 1 . / (4 . *BAL*SDK) ) 

BK-AK*BAR 

AKS«AK*AK 

BKS-BK*BK 

LIM-NSW+1 

IF (LIM . GT . 2 DM) LIM-IDM 
NPH-6 

NPH-2* (NPH/2) 

NPP-NPH+1 
PHA-. 0174533*20. 

DPH-PHA/NPH 

PH-.O 

DO  90  LPH-1,2 
WST-DPH*OST/ (3 . *PI ) 

SGN—  1. 

DO  80  IPH-1,NPP 

WF-3.+SGN 

IFdPH.EO.DMF-l. 

IF (IPH.EQ.NPP)WF-l . 

CPH-COS (PH) 

IF (IPH .GT . 1) GO  TO  40 
IF (LPH. GT. 1) GO  TO  40 
CPH-COS (DPH/10 . ) 

40  RS1-2.*AKS* (l.-CPH) 

RS2— AKS+BKS-2 . *AK*BK*CPH 
RH1-DSQRT (RSI) 

RH2-DSQRT (RS2 ) 

CALL  CIS I (CA, CIN, SA, RH1 ) 

CALL  C1SI (CB,CIN,SB,RH2) 

GI (1 ) — 2 . *CMPLX (CB-CA, SA-SB) 

DO  50  1-2, LIM 
DZ-DK* (1-1) 

DZS— DZ*DZ 
RA-DSQRT (RS1+DZS) 

RB-DSQRT (RS2+DZS) 

CALL  CISI (C1,CIN,S1,RA+DZ) 

CALL  CISI (C2, CIN, S2,RB+DZ) 

GP-CMPLX (C2-C1,  S1-S2) 

RAM— RS 1 / (RA+DZ ) 

RBM-R52/ (RB+DZ) 

CALL  CISI  (Cl, CIN,  SI, RAM) 

CALL  CISI (C2, CIN, S2.RBM) 

GM-CMPLX (C2-C1,  S1-S2) 

EGZ— CMP  LX  (COS  (DZ) ,  SIN  (DZ) ) 

50  GI (I) -GP*EGZ+GM/ECZ 

VJ(1)-VJ(1)+NF*NST*  (GI  (2)  -COK*GI  (I) ) 
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FRILLS. 2 


60 

78 

80 

90 


C 


IF (NEQ . LE . 1 ) GO  TO  78 
Kl-0 

DO  60  1-2, NEQ 
Kl-Kl+1 
K2-K1+1 
K3-K2+1 

IF  <K3 . GT . IDM) GO  TO  60 

GP-GI  (Kl) -2 . *CDK*GI (K 2)  +GI (K3) 

VJ  { I )  -VJ  (1 )  +WF*WST*GP 

CONTINUE 

SGN--SGN 

PH-PH4DPH 

DPH- <PI-PHA) /NPH 

PH-PHA 

VJ (1 ) -2 . *VJ (1) 

RETURN 

END 


7 


D-9 


■■ 


uuuuuuuuuuu 


c  TPU.l 

SUBROUTINE  TFLZ (C,U,VJ,W, Z, IER, IWR, 112, NEQ) 

MODIFIED  VERSION  OF  SUBROUTINE  FURNISHED  BY  CHARLES  KLEIN. 

SOLVES  SIMULTANEOUS  LINEAR  EQUATIONS. 

SET  IWR  -{1  OR  0)  TO  GET  (PRINTOUT  OR  NO-PRINTOUT). 

SET  112  -  1  ON  FIRST  CALL,  WHERE  MATRIX  INVERSION  IS  REQUIRED. 

SET  112  -  2  IF  MATRIX  Z  HAS  ALREADY  BEEN  INVERTED  ON  PREVIOUS  CALL. 
NEQ  -  NUMBER  OF  SIMULTANEOUS  LINEAR  EQUATIONS. 

Z(J>  IS  THE  FIRST  ROW  OF  THE  TOEPLITZ  IMPEDANCE  MATRIX 
VJ (J)  -  INPUT  VOLTAGE  COLUMN 
C(J)  -  OUTPUT  CURRENT  COLUMN 
U(J)  AND  W(J)  ARE  WORK  ARRAYS  OF  LENGTH  NEQ 
IF  IER  -  0  ,  NO  ERROR  OCCURRED 

COMPLEX  C()),D(1)  ,  VJ  (1),W(1),Z(1) 

COMPLEX  ALMDA, ALPHA, Cl , C2 , COEF , FAC, TAU1 ,  V, VI . V2 
2  FORMAT (1X,I5,F10.3,F15.7,F10.1) 

5  FORMAT  (1H0) 

IF(NEQ.GT.l)GO  TO  8 
C(1)-VJ(1) /Z<1) 

CNOR-CABS(CU)) 

GO  TO  100 

8  IF (112. NE . 1 ) GO  TO  45 

N-NEQ-1 
IER-0 

C  NORMALIZE  INPUT  MATRIX 
TAUl-Z(l) 

DO  10  II-1,N 
10  Z(II)-Z(I1+1) /TAU1 

ALMDA-1.-Z(1)*Z(1) 

U(l)»- 2(1) 

1-2 

15  KK-I-1 

ALPHA- ( . 0, . 0) 

DO  20  M-l.KK 
LL-I-M 

20  alpha-alpha+u(M)*z(ld 

ALPHA—  (ALPHA+2  (I)  ) 

ir (CABS (ALPHA) -EQ. . 0) GO  TO  130 

COEF -ALP HA /ALMDA 

ALMDA- ALMDA - COEF * ALPHA 

DO  30  «J—  1 ,  KK 

L-I-0 

30  W(J)-U(J)+COEF*U(L) 

DO  40  O-l, KK 
40  U(J)-W(J) 

U ( I ) -COEF 
IF (I .GE.N1GO  TO  45 
I-I  +  l 
GO  TO  15 

C  THE  FOLLOWING  COMPUTES  THE  ELEMENTS  OF  THE  INVERSE 
45  NH-(NEQ+l)/2 
FAC-ALMDA*TAU1 
NP-NEQ+ 1 
CNOR-  .  0 
DO  90  I-l.NH 
IF(I.NE.l)GO  TO  55 
W(l)  -1  .  /FAC 
DO  50  J-2.HEQ 
50  W(J)-U(J-1) /FAC 

GO  TO  TO 
55  Cl-U(I-l) 

NPI-NP-I 
C2-U (NPI ) 

DO  60  JJ— 1,N 

J-NP-JJ 

HPO-NP-J 
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60  W ( J) -W ( J-l ) + (Cl *0 ( J-l ) -C2*U (NPJ) ) /FAC 
Wfl)-U(I-l)  /FAC 
C  MATRIX  MULTIPLY 

to  v-(.o, .0) 
vi-(.o,  .0) 

DO  80  J-l, ME Q 
V2-VJ1J) 

V-V+V2-W  (J) 

NPJ-NP-J 

60  V1-V1+V2*W(NPJ) 

CU)-V 
NPI-NP-I 
C (NPI ) —VI 

IF  < IWR . LE . 0) GO  TO  90 
CA-CABS  (V) 

IF (CA . GT . CNOR) CNOR-CA 
CA-CABS (VI ) 

IF (CA . GT . CNOR) CNOR-CA 
90  CONTINUE 

100  IF (IWR. LE -  0) GO  TO  120 
C  PRINT  OUT  THE  SOLUTION  FOR  THE  CURRENTS  C ( J) 

WRITE  (6,  5) 

IF (CNOR . LE . 0 . ) CNOR-1 . 

DO  110  I-l.NEQ 
V-C(I) 

CA-CABS (V) 

CN-CA/CNOR 

PH-.O 

IF  (CA .  GT  .  0  . )  PH-51 .  29578*ATAN2  (AIMAG  (V)  ,  REAL  (V)  ) 
110  WRITE(6.2)I,CN,CA,PH 
WRITE (6, 5) 

120  RETURN 
130  IER-I 
RETURN 
END 
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TSPAR . 1 


SUBROUTINE  TSPAR!AK,DK,NEO,Z> 

TSPAR  sets  up  the  iapedance  aatrix  2<JI. 

Z(J)  -  First  tow  of  iapedance  aatrix  for  perfectly  conducting 
thin-wire  dipole  in  free  (pace,  using  Galerkin's  aethod  with 
overlapping  sinusoidal  basis  functions  and  aatching  the 
boundary  conditions  on  the  surface  of  the  wire. 

AK  -  k  *a ,  where  k  -  2'pl/laabde  and  a  -  wire  radius. 

DK  -  k'd,  where  d  -  segaent  length. 

NEO  -  nuaber  of  siaultaneous  linear  equations. 

REAL'S  DIS.RS 

COMPLEX  EID ( 90 ) , EM ( 90),EP!90),2(1) 

COMPLEX  CEM,CEP,END,EPD,EMD2,EPD2,Z11 ,222,G11 ,011 
DIMENSION  CID(90>,SID(90),CM(90),CP<90),SM(90),SP(90) 

DATA  GAM, ?2/. 51721 5664, 1.51019632/ 

DATA  ETA, PI/376. 727, 3. 14159/ 

IDM-90 

1  FORMAT! 3X, 'MUST  INCREASE  DIMENSIONS  IN  SUBROUTINE  TSPAR') 

2  FORMAT! 3X, 'ACTUAL  DIMENSION  I  DM  -  ',15, 6X, 

2  *  REQUIRED  DIMENSION  MAX2  -  ',15) 

I F! NEQ. LE.O) RETURN 
MAX2-NE0'2 
DO  19  I -1 , NEO 
14  Z<I)-!.0,.0) 

I F  t  MAX  2 . L  E . I  DM ) GO  TO  16 
WRITE! 6, 1 ) 

WRITE (6,2) 1DM , MAX2 
RETURN 

16  TDK-2. 'DK 
Sll-.O 
S13-TDK 
S2 1-DK 
S23-3 . *DK 
DO  20  N-1.MAX2 
I-N-l 
D2-1 'DK 
C I D ( N ) -COS ( OZ ) 

SIDIN)-SIN(DZ) 

20  E1D(N>-CMPLXICID(N),S1D(N)) 

CDK-COS ! DK ) 

SDK-SIN(DK) 

EPD— CMPLX ! CDK , SDK ) 

EMD-CMPLX ( CDK , -SDK  I 
EPD 2— EPD* EPD 
EMD2-EMD*EMD 
CEM-2.'CDK*EMD 
CEP-2 . 'CDK'EPD 
AK2-AK  *AK 

CSS-ETA/! 8 . 'PI'SOK'SDK ) 

NPH-6 

NPH-2' (NPH/2 I 

NPP-NPH'l 

PHA-. 0174533*20. 

DPH-PHA/NPH 

PH-.O 

DO  100  JPH-1,2 

CST-DPH'ETA/I 24 . *PI 'PI 'SDK'SDK ) 

C22-DPH/! 3 . 'PI ) 

SGN—  1 . 

DO  80  IPH-l.NPP 
CPH-COS ! PH ) 

SPH-S1N ! PH ) 

IP! IPH.GT.l (GO  TO  30 
IP! JPH.GT.l )GO  TO  30 
PH0-DPH/10. 

CPH-COS! PH0) 
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SPH-SIN(PHO) 

30  RH«AK*SPH 

RS-2.*AK2*(1.-CPHI 

RX-LSQRTIRS! 

Nr-3.«SGN 
irUPH.EO.nwr-1. 
i r  c iPH.EO.NPPjwr-i . 

WST-wr*CST 

W22-wr*C22 

DO  40  N-1.HAX2 

I-N-l 

DZ-I ‘DH 

DZS-DZ«DZ 

R-DSQRT(RS-DZS) 

ARG-R-OZ 

1MN.EQ.1  lARG-RX 

CALL  CISI(CP(N! ,CIN,SP(N> ,ARG) 

EP(N)-C«PLX(CP(N) ,-SP(N) ) 

IHN.GT.l IGO  TO  38 
CH(1)-CP(1) 

SM( 1 |-SP( 1 ) 

EM(l)-EPU) 

GO  TO  40 
36  ARG-RS/AFC. 

CALL  CISI (CH(N) ,CIN,S*fN) ,ARG) 

EM<N)-CHPLX(CH(N|.-SN(N) ) 

40  CONTINUE 

R-4.«(-CN(2)*2.*CP(l )-CP( 2) ) 
A*2.*CID(3J*(«CH(3)-2.*CH(2)*2.«CP(1)-2.*CH2)*CP(3)  ) 
B-2.*SID13)*(-SH(3)+2.*SH(2)-2.*SP(2)«SP(3> ) 
X-4.*(SH(2)-2.*SP(1)*SPI2) ) 
C-2.*CID(3)*(-SH(3)-2.*SH!2)-2.*SP(1 )-2 . ‘SP ( 2 ) -SP ( 3 > ) 
0*2.*SID(3)«(-CN(3)*2.*CR(2)-2.*CPI2)«CP(3) ) 
Z(1)-Z(1)*WST*CNPLX(R,X) 
ir<NEO.EO-l (GO  TO  70 

A-2.*CID<2>M-CH<3I  +  3.*CH<2>-4.*CK1)*3.*CP(2)-CP<3>) 
E-2.»SID(2I*(«SH(3)-2.*SN(2>*2.*SP(2)-SPI3>) 
P*CID(4)*(*CH(4I-2.«CH(3)*CN(2»*CP(2)-2.*CP(3)*CP(4)) 
G-»5ID(4>*(-SN(4>*2.*SMf  3)-SHf2)4SP|2)-2.*SP(3)»SP(4)  ) 
R-2.*CID(2)-(SN(3)-3.*SH(2)*4.*SP(1)-3.*SH2)«SP(3M 
H  «2. *SID(2)*(CH(3)-2.*C«(2)*2.*CP(2)-CP(3) ) 
I*CID(4)«(-SH(4)«2.*SN(3)-Sn(2)-SP(2)*2.*SP(3)-SP(4) ) 
J-»SID(4J*(-CB(4)*2.*CH(3>-CH(2>*CP12)-2.*CP(3)*CP(4>  > 
Z(2)-Z(2)«WST*CNPLX(R,X) 
ir(NEO.E0.2)GO  TO  70 
Sl-DK 

DO  60  N-3.NEQ 

Nl-N-1 

N2-N-2 

Nl-N-1 

N2-N-2 

CPA-CPt  N2 1-2.* CP (HI )«CP(N) 

CPB-2. *CP(N)-CP(N1 1-CPIN1 ) 

CPC-CP(N2)-2.*CP(N1 )*CP(N) 

CHA-C«<H2>-2. •CNIH1 )*CH(N) 

chb-2.*cimn)-ch<n1)-ch(M) 

CHC-CNIN2  >-2 . »CN(Nl )*CH(N) 

SPA-SP(H2)-2.*SP(Rl)+SP(N) 

SPB«2.*SP(N)-SP(H1 )-SP(Nl ) 

8PC-SP(N2)-2.«SP(Nl|*SP(N) 

SNA-SN (H2 )-2 . *5M( Ml )«SH( N) 

SHB-2.*SR(N)-SN(N1 ) -INI  Ml) 

SNC-SN ( N2 >-2 . •*« I  HI > *sn | N ) 

R-C2D(n2!*!CPA-C«»)*2.»ClD(N)«(CPB*CnB)-2.*SID(N)*tSPB-SHB) 
R  *CID(N2 )*ICPC»C«C)*SID(N2)»lSPC-SnCl 
iriN.CT. 3)r-r*sid(N2)*(spa-sra) 


u 
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TSPAR . 


TSPAR . 3 


X~CTO<N2>*iSPA+SNA)-2.*ClD(N)*(SPB*SNB)*2.*SlD(N)*(CPB-CKB) 
L  -CID(N2)MSPC*SHC)  +  SID(N2)MCPC-CHC> 

IF(N.GT.3>X-X  +  SID(H2)MCPA-CHA) 

60  X(N)«Z(N)*WST*CNPLX< R,X) 

70  PH-PH+DPH 
SO  SGN--SGN 

DPH-(P1-PBA)/NPH 
100  PH-PHA 
RETURN 
END 
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zsurf.i 


subroutine  zsuRriAX.crm.rnc.zs) 

ZSUrf  calculates  the  aurface  iapedance  IS  for  thin  wire. 
AX  -  k*a  where  k  -  2‘pi/laabda  and  a  ■  wire  radius, 
cun  •  conductivity  of  wire  (aegaahos/aeter ) 

CUM  -  -1  for  perfect  conductivity. 

EMC  -  frequency  (aegahertz). 

COMPLEX  BES.BES1.Z5 

DATA  ETA, SOT, TP/376. 727  27, 1.41421 356,6. 2831853/ 

SQSWE-1 . E6  *SQRT< CMN/TP/rnC/8 .85433  ) 

X- AR*SQSWE 
IP(X.GT.a. ) GO  TO  50 
T-X/8 . 

T2-T*T 

T4-T2*T2 

BER-<  < ( ( < (-.901E-5*T4-* . 122S52E-2 ) *T4- . 08349609 >*T4 
2*2. 641914  )*T4-32. 363456  >*T4  +  1 13. 7777  8 >*T4-6 4. )»T4«1. 

8EJ-< (III < ,11346£-3*T4-. 01103667 )*T4+. 5218561 5) *T4 
2-10.567658 )*T4+72. 817777 )*T4-11 3. 77778 )*T4*16. ) *T2 
BERP-X*T2*( (  (  ( <<-.394E-S*T4«. 4 5957 E- 3) *T4-. 02609253 )*T4 
2*. 66047849 ) *T4-6. 0681 481 )*T4* 14. 222222 )*T4-4. ) 

BEIP-XM ( 1 ( ( ( . 4609E-4*T4-. 3793B6E-2 )*T4».l 4677204 )*T4 
2-2.31167S1 )*T4*11. 377778 )*T4-10. 666667 )*T4*. 5) 
BES-CHPLXiBER,BEl ) 

BES1-.707107*CHPLX(BERP-BEIP,BERP«BEIP) 

GO  TO  100 

50  XP-.70710681*X 

XI- 1 ./X 

r-( <-.04  59205*Xl*.390625E-2)*Xl*.088  388  35)*XUl . 

T-< (-. 04603559*X1-. 0625 )*X1-. 068388  35  I *X1- . 39269907*XP 
BES-r*CNPLX(COS(T) ,SIN(T) ) 

r-( ( . 11290231*X1*. 0351 562 5 )‘X1-. 26516505 ) *X1* 1 . 

T-( ( .1360097*Xle,1875)*Xl* . 2651 6505 ) *X1*1 . 1780972+XP 
BES1-P •CMPLXI COS ( T ) , SIN( T } ) 

100  ZS— CMPLX(  1 .  .-1 .  )  *ETA*BES/BES1/SQT/SQSWE 
RETURN 
END 


Appendix  II.  Output  data  from  flsit-Surteea .For 


AL 

0.00300 


BL 

0.10000 

0.20000 

0.30000 

0.40000 

0.50000 

0.60000 

0.70000 

0.80000 

0.90000 

1.00000 

1.10000 

1.20000 

1.30000 

1.40000 

1.50000 

1.60000 


HL 

0.22900 


RFIN 

43.93178 

28.43193 

25.36656 

29.26224 

34.19707 

35.41125 

32.85429 

30.25837 

30.52595 

32.74195 

33.91415 

32.81260 

31.15025 

31.06938 

32.45232 

33.36460 


RINF 

32.22021 


XFIN 

-27.96955 

-20.13158 

-11.10953 

-5.85667 

-6.53444 

-10.37363 

-12.58671 

-11.25765 

-8.70713 

-8.13379 

-9.83655 

-11.36605 

-10.87897 

-9.27844 

-8.72108 

-9.77519 


14* 


XINF 

-9.93279 
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